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Abstract 

We analyze rigorously error estimates and compare numerically spatial/temporal resolution of various nu¬ 
merical methods for the discretization of the Dirac equation in the nonrelativistic limit regime, involving 
a small dimensionless parameter 0 < e <C 1 which is inversely proportional to the speed of light. In this 
limit regime, the solution is highly oscillatory in time, i.e. there are propagating waves with wavelength 
0(s 2 ) and 0(1) in time and space, respectively. We begin with several frequently used finite difference time 
domain (FDTD) methods and obtain rigorously their error estimates in the nonrelativistic limit regime by 
paying particular attention to how error bounds depend explicitly on mesh size h and time step r as well 
as the small parameter e. Based on the error bounds, in order to obtain ‘correct’ numerical solutions in 
the nonrelativistic limit regime, i.e. 0 < £ C 1, the FDTD methods share the same e-scalability on time 
step and mesh size as: r = 0(e 3 ) and h = 0{y/e). Then we propose and analyze two numerical methods 
for the discretization of the Dirac equation by using the Fourier spectral discretization for spatial deriva¬ 
tives combined with the exponential wave integrator and time-splitting technique for temporal derivatives, 
respectively. Rigorous error bounds for the two numerical methods show that their e-scalability is improved 
to r = 0(e 2 ) and h = 0(1) when 0 < e <C 1. Extensive numerical results are reported to support our error 
estimates. 

Keywords: Dirac equation, nonrelativistic limit regime, finite difference time domain method, exponential 
wave integrator spectral method, time splitting spectral method, e-scalability 


1. Introduction 

The Dirac equation, which plays an important role in particle physics, is a relativistic wave equation 
derived by the British physicist Paul Dirac in 1928 ESI Eg imug. It provided a description of elementary 
spin-1/2 massive particles, such as electrons and positrons, consistent with both the principle of quantum 
mechanics and the theory of special relativity. It was the first theory to fully account for relativity in the 
context of quantum mechanics. It addressed the fine details of the hydrogen spectrum in a completely 
rigorous way and predicted the existence of a new form of matter, antimatter (4|. Since the graphene 
was first produced in the lab in 2003 [U EH [621 IMl US], the Dirac equation has been extensively adopted 
to study theoretically the structures and/or dynamical properties of graphene and graphite as well as two 
dimensional (2D) materials [60 j . This experimental advance renewed extensively the research interests on the 
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mathematical analysis and numerical simulations of the Dirac equation and/or the (nonlinear) Schrodinger 
equation without/with external potentials, especially the honeycomb lattice potential [3, [35]. 

We consider the three dimensional (3D) Dirac equation for describing the time evolution of spin-1/2 
massive particles, such as electrons and positrons, within external time-dependent electromagnetic potentials 

Ha m 


ihd t 'f>(t,x) = —ich'y^ajdj 
j=i 


d'(t,x) + e V(t, x )/4 — ^ Aj(t,x.)aj \I/(f,x). (1.1) 


i= i 


Here, i = \J— 1, t is time, x = (aq, x 2 , x^) T £ M 3 (equivalently written as x = {x,y,z) T ) is the spatial 
coordinate vector, dk = {k = 1,2,3), := 4 /(£,x) = x), i/> 2 (t, x), ^(t, x), ip 4 (t, x)) T £ C 4 is 

the complex-valued vector wave function of the “spinorfield”. I n is the n x n identity matrix for n £ N, 
V := V(t,x) is the real-valued electrical potential and A := A(£, x) = x), A 2 {t 1 x), As(t, x)) T is the 

real-valued magnetic potential vector, and hence the electric field is given by E(i, x) = — W — d t A and the 
magnetic Held is given by B(t, x) = curl A = V x A. The physical constants are: c for the speed of light, 
m for the particle’s rest mass, H for the Planck constant and e for the unit charge. In addition, the 4x4 
matrices aq, a 2 , a. 3 and /3 are defined as 


Oil = 


0 G 1 
G 1 0 


a 2 = 


0 a 2 

G-2 0 


C*3 = 


0 tx 3 
cr 3 0 


p = 


I 2 

0 


0 

—I2 


( 1 . 2 ) 


with g !, g 2 , CT 3 (equivalently written g x1 g v , g z ) being the Pauli matrices defined as 


fTl = 


0 1 
1 0 


02 = 


0 -* 
i 0 


03 


1 0 
0 -1 


(1.3) 


In order to scale the Dirac equation (1.1 1 , we introduce 


t = ^~, ^(£,x) = x 3 J 2 T(f,x), V (t, x) = V ^’ , Aj (t, x) = 

Tc (Z/o /1c si a 


3 = 1,2,3, (1.4) 


where x a , £ s and A s are the dimensionless length unit, time unit and pote ntial unit, respectively, satisfying 

t„x 3 , /2 


and A s = rn ^- with v = ^ being the wave speed. Plugging (1.41 into (1.1 


t s = 

and then removing all” we obtain the following dimensionless Dirac equation in' 


, multiplying by 


3D 


h ’ 


. 3 1 3 

id t ^(t,x)= - - ^ ajdj + '&(t, x) + H(£,x)I 4 - Aj(t,x)atj ^(£,x), x£ 

£ i =1 £ 3 =1 


where e is a dimensionless parameter inversely proportional to the speed of light given by 


an v 

0 < e := —— = — < 1. 
t„ c c 


(1.5) 


( 1 . 6 ) 


We remark here that if one chooses the dimensionless length unit x s = — , t s = — and A s = m£ - 


in (h.4|, then e = 1 in (|1.6|) and Eq. (1.5) with e = 1 takes the form often appearing in the literature 


El 11711211 [23 ; 33. [391 1471 151] . This choice of x s is appropriate when the wave speed is at the same order 
of the speed of light. However, when the wave speed is much smaller than the speed of light, a different 
choice of x s is more appropriate. Note that the choice of x s determines the observation scale of the time 
evolution of the particles and decides: (i) which phenomena are ‘visible’ by asymptotic analysis, and (ii) 
which phenomena can be resolved by discretization by specified spatial/temporal grids. In fact, there are 
two important parameter regimes: One is e = 1 (<t=> x s = t s = and A s = 2 ^-), then Eq. (1.51 


describes the case that wave speed is at the same order of the speed of light; the other one is 0 < e <C 1 , 


then Eq. (1.5) is in the nonrelativistic limit regime. 
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Similarly to the dimension reduction of the nonlinear Sclirodinger equation and/or the Schrodinger- 
Poisson equations with/witlrout anisotropic external potentials |8j, when the initial data ^((hx) and the 
electromagnetic potentials V ( t , x) and A(f, x) are independent of z and thus the wave function ’I' is formally 
assumed to be independent of z, or when the electromagnetic potentials P(f,x) and A(f,x) are strongly 
confined in the ^-dir ection and thus T is formally assumed to be concentrated on the xy-plane, then the 3D 
Dirac equation (1.5) can be reduced to the Dirac equation in 2D with x = ( x,y) T G M 2 as 


id t ^(t,x) = + -^0 ^{t,x) + V ( t , x) J 4 - ^ Aj ( t , x)ctj 4/(f,x), xG 

£ j =i £ j =i 


(1.7) 


This 2D Dirac equation has been widely used to model the electron structure and/or dynamical properties 
of graphene since they share the same dispersion relation on the Dirac points [B [HU EH [Ml ESI EH EH EZ] - 
Similarly, under the proper assumptions on the initial data and the external electromagnetic potential, the 
3D Dirac equation (1.51 can be reduced to the Dirac equation in ID with T = ^(f, x) as 


id t $>(t,x) = 


i q , 1 
— ot\d x H—,, 


T(t,x)+ V{t 7 x)I± — Ai(t,x)ai 4 /(t,x), itl. 


( 1 . 8 ) 


In fact, the Dirac equation in 3D (1.5), in 2D (1.7) and in ID (1.8) can be written in a unified way in 
d-dimensions (d = 1, 2, 3) 


id t ^(t,x) = --"Y^ajdj + ^/3 ^(t,x) + P(f, x)/ 4 - ^ A, (f, x)ay 4>(f,x), xG 
£ j =l £ i=i 


and the initial condition for dynamics is given as 

^(f = 0, x) = T 0 (x), 


x G 


(1.9) 


( 1 . 10 ) 


The Dirac equation (1.9) is dispersive and time symmetric. Introducing the position density pj for the 
j-component (j = 1,2, 3,4) and the total density p as well as the current density J(f,x) = (Ji(f, x), J 2 {t,x ), 


p(t,x) = ^ft(t,x) = Pj(t,x) = |i/’j(t,x)| 2 , l<j<4; Ji(t,x) = 1 = 1,2,3, (1.11) 

i=i e 

— —T 

where / denotes the complex conjugate of / and ^ , then the following conservation law can be 

( 1 . 12 ) 


obtained from the Dirac equation (1.9) 


d t p(t, x) + V • J(f, x) = 0, x G R , t > 0. 


Thus the Dirac equation (1.9) conserves the total mass as 

4 


ll^(^-)|| 2 := / \'S>{t,x)\ 2 dx= f J2\ipj{t,x)\ 2 dx= ||T(0,-)|| 2 = H’holl 2 , t > 0. (1.13) 

JR d JR d . =1 


If the electric potential V is perturbed by a real constant V°, e.g. V(t,x) —> V(t,x) + V 0 , then the 
solution 4/(f,x) —>■ e~ lV t 'Sf(t,x) which implies the density of each component pj (j = 1,2,3,4) and the 
total density p unchanged. When d = 1, if the magnetic potential A\ is perturbed by a real constant A®, 
e.g. Ai(t,x) —> Ai(t,x) + A°, then the solution 4>(t,x) —> e l2lltQl 4'(t,x) which implies the total density p 
unchanged; but this property is not valid when d = 2, 3. In addition, when the electromagnetic potentials are 
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time-independent, i.e. V(t, x) = V(x) and Aj(t. x) = Aj (x) for j = 1,2,3, the following energy functional 
is also conserved 


d d 

E(t):= I ( + ] dx=E(0), t > 0. (1.14) 

' l=i £ l=i 


lR d 


Furthermore, if the external electromagnetic potentials are const ants , i.e. V(t,x) = V° and A,(f,x) = 
Aj for j = 1,2,3 with A° = (A) 1 ,..., A°) T , the Dirac equation (1.91 admits the plane wave solution as 
^(f, x) = g e i( k x -w*), where the time frequency w, amplitude vector B G R 4 and spatial wave number 
k = (ki, ..., kd) T G satisfy the following dispersion relation 


wB = 


a 

[E 

l=i 


‘-f- A °)^ + ^P + V Ul 4 


B, 


(1.15) 


which immediately implies the dispersion relation of the Dirac equation (1.9) as 


:= w(k) = V° ± \\! 1 + £ 2 |k — eA 0 !" 2 , k G 


(1.16) 


Plugging (1.2) and (1.3) into (1.7), the 2D Dirac equation (1.7) can be decoupled as 

i 1 

idt^i = --(d x -id y )ip 4 + + V(t,x)il >i - [Ai(t,x) - iA 2 (t,x)\ 0 4 , 

X 1 

id t 4 >4 = --{d x +id y )ip i- + P(t,x)0 4 - [Ai(t,x) + iA 2 (t,x)\ 0i, x G 

£ 

X 1 

id t if >2 = --{d x +id v )ip 3 + ^0 2 + V(i,x)0 2 - [Ai(t,x) + L4 2 (f,x)] 0 3 , 

X 1 

idti >3 = -- (0® - idy)ip 2 -y03 + V(i,x)0 3 - [Ai(t,x) - L4 2 (f,x)] 0 2 , x G 


(1.17) 


(1.18) 


Eq. (1.18 ) wil l collapse to ( 1.17| ) under the transformation y —1 — y and A 2 —1 — A 2 . Thus, in 2D, the Dirac 
equation (1.7) can be reduced to the following simplified PDEs with 4> = <b(f,x) = (</>i(f, x), </> 2 (f,x)) T G C 2 

id t $(t,x) = (cti<9 x + <j 2 d y ) + -^cr 3 ]<F(t,x) + \v(t,x)I 2 - A 1 (t,x)o 1 - A 2 (t,x)a 2 j$(t,x), x G R 2 , 

L e e z J L J 

( 1 -!9 ) 

where $ = (0i,0 4 ) T (or $ = ('02,'*/’ 3 ) T under the transformation ?/ —> —y and A 2 —> —A 2 ). Similarly, 
in ID, the Dirac equation ( |1.8| ) can be reduced to the following simplified PDEs with $ = = 

(cj) 1 (t,x),(t> 2 {t,x)) T 


id t ${t,x) = 


-~oid x + —.03 


$(t,x)+ V(t,x)I 2 — Ai(t,x)ai i£l, (1-20) 


where $ = (0i,04) T (or $ = (02,'*/’ 3 ) T ). Again, the Dirac equation in 2D (1.19) and in ID (1.20) can be 
written in a unified way in d-dimensions (d = 1, 2) 


id t $(t,x) = -- Vjdj + 9 o ~3 x) + V(t,x)I 2 -'^Aj(t,x)<j j $(t,x), x € K d , (1.21) 
£ i=i £ 3 = 1 


and the initial condition for dynamics is given as 

<h(f = 0,x) = $ 0 (x), x G R d . 


( 1 . 22 ) 
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The Dirac equation (1.211 is dispersive and time symmetric. By introducing the position density pj for the 
j-th component (j = 1, 2) and the total density p as well as the current density J(t, x) = (Ji(t, x), J 2 (t, x)) T 


2 1 
p{t,x) = ^Pj&x) = Pi(*) x ) = l^( t > x )| 2 » Jj(i,x) =i = 1,2, 


(1.23) 


1=1 


the conservation law (1.12) is also satisfied [23]. In addition, the Dirac equation (1.21) conserves the total 
mass as 


ll$(V)|| 2 := / |$(i,x)| 2 dx= f ^|<() j (t,x)| 2 dx= ||$(0,-)f = ll^of, t > 0. (1.24) 

JR d JR d - =1 

Again, if the electric potential 14 is perturbed by a real constant 14°, e.g. 14(t,x) -> V(t,x) + V°, the 
solution <F(£,x) —► e _lV *<!>(£, x) which implies the density of each component pj ( j = 1,2) and the total 
density p unchanged. When d = 1, if the magnetic potential A\ is perturbed by a real constant A®, e.g. 
Ai(t,x) —► Ai(f,x) + Ai, the solution $(£, x) —► e lA A ai $(£,x) implying the total density p unchanged; 
but this property is not valid when d = 2. When the electromagnetic potentials are time-independent, i.e. 
14(f,x) = 14(x) and A ;i (t. x) = A.y(x) for j = 1,2, the following energy functional is also conserved 


E{t) := [ (--V$V i 9 i $+4 $ * ct 3$ + 1 / (x)|$| 2 -y^A,-(x)$Vj$ ) dx=£(0), t > 0. (1.25) 

■ / * i V £ i=! £ J =1 / 

Furthermore, if the external electromagnetic potentials are constants, i.e. V(f, x) = V° and Aj(t,x ) = A° 
for j = 1,2, the Dirac equation (1.211 admits the plane wave solution as $(t,x) = B e i(k.x-a,i) j where the 


time frequency u>, amplitude vector B S 
following dispersion relation 


and spatial wave number k = (fci,..., kd) T G satisfy the 


jB = 


d , 

£ 

l=i v 


£ J 


4-+ V°I 2 


B. 


(1.26) 


which again implies the dispersion relation ( 1.16| ) of the Dirac equation (1.21) for rf = 2,1. 

For the Dirac equation (1.9) with e = 1, i.e. 0(l)-speed of light regime, there are extensive analytical 
and numerical results in the literatures. For the existence and multiplicity of bound states and/or standing 
wave solutions, we refer to (31ES M EU EH [74j and references therein. For the analysis of the clas- 
sical/semiclassical limits via the Wigner transform techniques, we refer to [H (7, HI [HI ESI 'tL S3 and 
references therein. For the numerical methods and comparison such as the finite difference time domain 
(FDTD) methods and the Gaussian beam methods, we refer to [51 [7,178,177,175, [27, [581 [754 H3] and references 
therein. However, for the Dirac equation (1.9) with 0 < e 1, i.e. nonrelativistic limit regime (or the scaled 
speed of light goes to infinity), the analysis and efficient computation of the Dirac equation (|1.9[) (or (1.21)) 


are mathematically rather complicated. The main difficulty is due to that the solution is highly oscillatory 
in time and the corresponding energy functionals (1.14) and (1.25) are indefinite in ei and become un¬ 
bounded when e —► 0. There are extensive mathematical analysis of the (semi)-nonrelativistic limit of the 


Dirac equation (1.9) to the Pauli equation [51 [H (T5J [H EH1IH El El El El El and/or the Schrodinger 
equation when e —> 0 [18]. These rigorous analytical results show that the solution propagates waves with 
wavelength 0(s 2 ) and 0(1) in time and space, respectively, when 0 < £ <C 1. In fact, the oscillatory structure 
of the solution of the Dirac equation (1.9) when 0 < £ <C 1 can be formally observed from its dispersion rela¬ 
tion (1.15) (or ( |1.26| )). To illustrate this further, Figure | lT| shows the solution of the Dirac equation (1.21) 

with d = 1, V(t,x) = yqE/&, A\{t,x) = and 4> 0 (a;) = (ex p(—z 2 /2),exp (—{x — 1) 2 /2)) T for different e. 

This highly oscillatory nature of the solution of (1.91 (or (1.21)) causes severe numerical burdens in practical 


computation, making the numerical approximation of (1.9) (or (1.21)) extremely challenging and costly in 
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Figure 1.1: The solution <j>i(t = \,x) and <t>\(t,x 
e. Re(/) denotes the real part of /. 


0) of the Dirac equation (1.211 with d = 1 for different 


the nonrelativistic regime 0 < £ C 1. In m , the resolution of the time-splitting Fourier pseudospectral 
(TSFP) method was studied for the Maxwell-Dirac equation in the nonrelativistic limit regime. 

Recently, different numerical methods were proposed and analyzed for the efficient computation of the 
Klein-Gordon equation in the nonrelativistic limit regime HU H21 [34] and/or highly oscillatory dispersive 
partial differential equations (PDEs) 0, 10, 13, H4] , To our knowledge, so far there are few results on the 
numerics of the Dirac equation in the nonrelativistic limit regime. The aim of this paper is to study the effi¬ 
ciency of the frequently used FDTD and TSFP methods applied to the Dirac equation in the nonrelativistic 
limit regime, to propose the exponential wave integrator Fourier pseudospectral (EWI-FP) method and to 
compare their resolution capacities in this regime. We start with the detailed analysis on the stability and 
convergence of several standard implicit/semi-implicit/explicit FDTD methods [70j. Here we pay particular 
attention to how the error bounds depend explicitly on the small parameter e in addition to the mesh size 
h and time step r. Based on the estimates, in order to obtain ‘correct’ numerical approximations when 
0 < e « 1, the meshing strategy requirement (e-scalability) for those frequently used FDTD methods is: 
h = O(^fe) and r = 0(e 3 ), which suggests that the standard FDTD methods are computationally expensive 
for the Dirac equation (1.9) as 0 < e <C 1. To relax the e-scalability, we then propose the EWI-FP method 
and compare it with the TSFP method, whose e-scalability are optimal for both time and space in view of 
the inherent oscillatory nature. The key ideas of the EWI-FP are: (i) to apply the Fourier pseudospectral 
discretization for spatial derivatives; and (ii) to adopt the exponential wave integrator (EWI) for integrat¬ 
ing the ordinary differential equations (ODEs) in phase space [401 fig] which was well demonstrated in the 
literatures that it has favorable properties compared to standard time integrators for oscillatory differential 
equations [40, 48} 09} [50]. Rigorous error estimates show that the e-scalability of the EWI-FP method is 
h = 0(1), and r = 0(e 2 ) for the Dirac equation with external electromagnetic potentials, meanwhile, the 
e-scalability of TSFP method is h = 0(1) and r = 0(e 2 ). Thus, the EWI-FP and TSFP offer compelling 
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advantages over commonly used FDTD methods in temporal and spatial resolution when 0 < e 1. 

The rest of this paper is organized as follows. In Section [2j several second-order FDTD methods are 
reviewed and their stabilities and convergence are analyzed in the nonrelativistic limit regime. In Section 
[3j an exponential wave integrator Fourier pseudospectral method is proposed and analyzed rigorously. In 
Section |4j a time-splitting Fourier pseudospectral method is reviewed and analyzed rigorously. In Section [5j 
numerical comparison results are reported. Finally, some concluding remarks are drawn in Section [hj The 
mathematical proofs of the error estimates are given in the appendices, where extensions of EWI-FP and 
TSFP to higher dimensions are also presented. Throughout the paper, we adopt the standard notations 
of Sobolev spaces, use the notation p < q to represent that there exists a generic constant C which is 
independent of h , r and e such that |p| < C q. 


2. FDTD methods and their analysis 


In this section, we apply the commonly used FDTD methods to the Dirac equation (1.9) (or (1.21)) and 
analyze their stabilities and convergence in the nonrelativistic limit regime. For simplicity of notations, we 


shall only present the numerical methods and their analysis for (1.21) in ID. Generalization to (1.9) and/or 
higher dimensions is straightforward and results remain valid without modifications. Similarly to most works 
in the literatures for the analysis and computation of the Dirac equation (cf. 17, 23 SSi 323 EH IMS E7] EZSJ 
and references therein), in practical computation, we truncate the whole space problem onto an interval 
Q = (a, b) with periodic boundary conditions, which is large enough such that the truncation error is 


negligible. In ID, the Dirac equation (1.21) with periodic boundary conditions collapses to 


id t ${t,x) = 


1 




~-aid x + —a 3 
I £ £ z 

<F(t,a) =$(t, b), d x $(t,a) = d x $(t,b), t > 0, 


V(t,x)I 2 — Ai(t,x)ai $(t,x), itll, t > 0, (2.1) 

$(0,*) = $ o (a0, xeU, (2.2) 


where <J> 0 (a) = $o(&) and $o(a) = $q(&). 


2.1. FDTD methods 

Choose mesh size h := Ax = with M being an even positive integer, time step r := At > 0 and 
denote the grid points and time steps as: 


Xj := a + jh , j = 0,1,..., M; t n := nr, n = 0,1,2 ,.... 

Denote Xm = {U = (Uq, U\,..., Um ) 1 \ Uj £ C 2 , j = 0,1,..., M, Uq = Um} and we always use U- 1 = Um -i 
and Um+i = U\ if they are involved. For any U £ Xm, we denote its Fourier representation as 


M/2-1 M/2-1 

Uj = Uie iMxj ~ a) = Uie 2ijlv/M , 


j = 0,1,..., M, 


l=-M /2 


1— — M/2 


where /q and Ui £ C 2 are defined as 


m = 


2 hr 
b — a’ 


M-l 


U = j- V Uj I = _ ^ M - L 

M J ’ o ’ ’ o 


3= 0 


The standard t 2 -norm in Xm is given as 


M-l 


\mi = hj2 \u,\ 2 , u£x m . 

3=0 


(2.3) 


(2.4) 


(2.5) 


Let be the numerical approximation of $(t n ,Xj) and VJ 1 = V(t n ,Xj), V ? ” +1 ^ 2 = V(t„ +r/2,Xj), A" • = 
Ai(t n , Xj) and A ^ 1 ^ 2 = Ai(t n + r/ 2, Xj ) for 0 < j < M and n > 0. Denote <h ra = ($q, • • •, r } iI ) T £ X M 
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as the solution vector at t = t n . Introduce the finite difference discretization operators for j = 0,1 ,M 
and n > 0 as: 


$"+' - <j>« 
§ +^u = - 3_ 

1 J -T 






2r 


= 


2/i 


®"+l i 
J ^ J 


(E> 2 = 

3 


Here we consider several frequently used FDTD methods to discretize the Dirac equation (2.1) for 

j = 0,1,..., M — 1. 

I. Leap-frog finite difference (LFFD) method 


= 


* r 1 
-Cl Ox d—yC3 




vpi 2 - 


T”, n > 1. 


( 2 . 6 ) 


II. Semi-implicit finite difference (SIFD1) method 

a i $”+ 1 $ n_1 

iS t $? = --a,S x <$] + -^o 3 ^ - - 

J e J e 2 - 2 


Vfh - A^ox 


$n+l + <j>n-1 

— --—-—, n> 1. 


(2.7) 


III. Another semi-implicit finite difference (SIFD2) method 



i 1 1 

a n+1 -j- <h n_1 


= 

-H-o ^3 

£ £ z 

3 3 | 

2 

[vrh-A-.o i] 


n > 1. 


( 2 . 8 ) 


IV. Crank-Nicolson finite difference (CNFD) method 


i6+<S> 


n 

3 


- (TlS x + 

L £ 



$ n+1 / 2 

3 


t rn+l/2 r ^n-f-1/2 

Fj t 2 - 


ci 


$ 


1+1/2 


n > 0. 


(2.9) 


The initial and boundary conditions in (2.2) are discretized as: 

p.n+1 


$M +1 = 


QU +1 = Qn+l^ 


^ = ^ 0 (Xj), j = 0, 1, ..., M. 


( 2 . 10 ) 


Using Taylor expansion and noticing (2.1), the first step for the LFFD (2.6), SIFD1 (2.7) and SIFD2 ( |2.8| ) 
can be computed as 


$j=^+T 


1 

r 


sin oi& 0 (xj) ~i{~ sin ) c 3 + V°I 2 - A^jOr 




j = 0 , 1 ,-.., M. (2.11) 


In the above, we adapt h sin (^) and i sin (p-) instead of j and ^ such that (2.11) is second order in term 


of r for any fixed 0 < e < 1 and ||$ 1 || oo := max |$)| < 1 for 0 < e < 1. We remark here that they can be 

0<j<M J 

simply replaced by 1 when e = 1. 

The above four methods are all time symmetric, i.e. they are unchanged under r ++ —r and n + 1 ++ n— 1 
in the LFFD, SIFD1 and SIFD2 methods or n + 1 ++ n in the CNFD method, and the memory cost is the 


same at O(M). The LFFD method (2.6) is explicit and its computational cost per step is O(M). In fact, 
it might be the simplest and most efficient discretization for the Dirac equation when e = 1 and thus it has 


been widely used when e = 1. The SIFD1 method (2.7) is implicit, however at each time step for n > 1, the 


corresponding linear system is decoupled and can be solved explicitly for j = 0,1,..., M — 1 as 


$"+ 1 

3 


(i - TV n )I 2 - ^<73 + TA'l 0 + 


((* + tU")/ 2 + ^2 c73 - tA ijOi^ 




( 2 . 12 ) 


and thus its computational cost per step is 0(M). 





































The SIFD2 method ( |2.8| ) is implicit, however at each time step for n > 1, the corresponding linear system 
is decoupled in phase (Fourier) space and can be solved explicitly in phase space for l = —M/2,..., M/2 — 1 
as 


,, , in i rsin(m/i) t 

= I ' h - -vsP" 1 - ^ 


Tsm(nih) 

ih H- 7 - CT i 

eh 



), + 2 r(G"$"), 


, (2.13) 


where G" = (Gg, G",..., G)^ ) T £ Am with G™ = — A^-ai + V ^ /2 for j = 0,1,..., M, and thus its 


computational cost per step is 0(M In M). The CNFD method (2.9) is implicit and at each time step for 


n > 0, the corresponding linear system is coupled and needs to be solved via either a direct solver or an 
iterative solver, and thus its computational cost per step depends on the linear system solver, which is 
usually much larger than O(M), especially in 2D and 3D. Based on the computational cost per time step, 
the LFFD method is the most efficient one and the CNFD method is the most expensive one. 


2.2. Linear stability analysis 

In order to carry out the linear stability analysis for the FDTD methods via the von Neumann method 
123 . we assume that A\{t,x) = A® and V{t, x) = V° with A® and V° being two real constants in the Dirac 
equation (2.1). Then we have the following results for the FDTD methods: 


Lemma 1. (i) The LFFD method (2.6) is stable under the stability condition 

e 2 h 


0 < t < 


| V°\e 2 h + \/h 2 + e 2 (l + e/i|y4°|) 2 ’ 

(ii) The SIFD1 method \2. 7| ) is stable under the stability condition 

0 < t < eh, h > 0, 0 < e < 1. 


h > 0, 0 < e < 1. 


(Hi) The SIFD2 method (2.8) is stable under the stability condition 

1 


0 < T < 


h > 0, 0 < £ < 1. 


(2.14) 


(2.15) 


(2.16) 


in + w 

(iv) The CNFD method {2. ,9[ j is unconditionally stable, i.e. it is stable for any r,h > 0 and 0 < £ < 1. 
Proof: (i) Plugging 


M/ 2-1 M/ 2-1 

^ ( 5 °); ^ £” ( 5 °); e 2i ^ ln / M , 

Z=-M/2 Z=-M/2 


j = 0,1,..., M, n > 0, 


(2.17) 


with £ C and ('ll 0 ); being the amplification factor and the Fourier coefficient at n = 0, respectively, of the 
Z-tfi mode in the phase space into (2.6), using the orthogonality of the Fourier series, we obtain 


{it - 1 )I 2 - 2iri ( A°a! - V°I 2 -^a 3 - 


sin (jm/i) 

eh 


-ay 


M M 

= 0, =-,...,-1. 

2 2 


Substituting (1.3) into (2.18), we get that the amplification factor i satisfies 

$-21x0^1-1 = 0, Z = -y,...,y- 1, 


(2.18) 


(2.19) 


where 


0i = -V Q ± h 2 + £ 2 (A^eh - sin {/Jih)) 2 , l 


M M 
~~^ 2 ''"'^ 2 ~ 
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Then the stability condition for the LFFD method (2.6 1 becomes 


M M 

161 <1 *=* Kl< 1, / = - 


( 2 . 20 ) 


which immediately implies the condition (2.14). 

(ii) Similarly to (i), plugging (2.17) into the SIFD1 method (|2.7|), we have 


( 6 2 - 1 )h + 1 ) A'ia, V"I 2 - + 


1 \ 2*r6 sin(/x ; /i) 


eh 


-o’ 1 


= 0, Z = (2.21) 


Noticing (1.3), under the condition (2.15), we can get |6| < 1 for l = — ,..., ^ — 1, and thus it is stable, 

(iii) Similarly to (i), plugging ( 2.1 7| ) into the SIFD2 method (2.8), we have 


( 6 2 — 1)^2 + * t ( 6 2 + 1 ) ( + 


sin (pih) 
eh 


cri - 2ir^ l (A° 1 a 1 - V°I 2 ) 


= 0, Z = -^..,^-l. (2.22) 


Noticing (1.3), under the condition (2.16), we obtain 


161 <i, l = 


and thus it is stable. 


(iv) Similarly to (i), plugging (2.17) into the CNFD method (2.9), we obtain 


(6 — 1)^2 + tt (6 + 1) ( ~2 a 3 ~ Ai<Ji — V°I-2 + 


si n (mh ) 
eh 


-(Jl 


„ , M M „ 

-°, l-~ —-i. 


Noticing (1.3), we have for l = —— 1, 


161 = 


2 + irdi 


1 

e 2 h 


= 1, 0i = V° ± -A^\/h 2 + e 2 {A\eh — sin (pih)Y 


2 — irOi 

Thus it is unconditionally stable. 

2.3. Mass and energy conservation 

For the CNFD method (2.9), we have the following conservative properties. 

Lemma 2. The CNFD \2.9\ ) conserves the mass in the discretized level, i.e. 

M -1 M—l M—l 

m\p ■■= h E ! $ ”| 2 = h E l^l 2 = H $ °IIp = h E l $ o(^-)| 2 , n > 0. 


(2.23) 

(2.24) 

□ 


(2.25) 


l=o 


1=o 


l=o 


Furthermore, ifV[t,x) = V(x) and Ai(t,x) = A\{x) are time independent, the CNFD (2.9) conserves the 
energy as well, 


., M—l , M—l M—l M—l 

K = - V E ($i) + -2 E (*")**3*? +E -hJ2 A lt j (§>*)* ai&j 


l=i 


l=o 


1=0 


1=0 


(2.26) 


=E°, n> 0 , 


where Vj = V(xj ) and .A^- = ^(ar,-) /or j = 0,1,..., M. 
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Proof: (i) Firstly, we prove the mass conservation (2.25). Multiplying both sides of (2.9) from left by 
hr ($" +1 / 2 )* and taking the imaginary part, we have 


hi $” + i i 2 = h \^\ 2 - ^ [(^+ i / 2 ) v 1 4^" +1/2 + ($; + 1 / 2 ) T CTl 4 ^ 


—n+i/ 2 ' 


j = 0,1,..., M — 1. (2.27) 


Summing (2.27) for j = 0,1,..., M — 1 and noticing (1.3), we get 

M -1 


|$«+ 1 ||2 2 = ||$"||2 


rh 

2 e 


£ [(*y ,/2 r *A*r /2 +(*r ,/ y ^**J +1/2 


1=0 

M—l 


= 11$ 


n || 2 


E [« 1/2 r ^;i 1/2 + (*r l/2 ) T 


2 e 


1=o 


-(C7T -i$r l/2 - WXiT 


= n> 0, 


(2.28) 


which immediately implies (2.25) by induction. 

(ii) Secondly, we prove the energy conservation (2.26). Multiplying both sides of (2.9) from left by 
2 h ($” +1 — $"■)* and taking the real part, we have 


— h Re 


-($? +1 - $, n )Vi^($; +i + $7) 


+4 m +l r^r - 


+ hV 3 {\^ +1 \ 2 - |$”| 2 ) - hA ltj [($” +1 )Vi$" +1 - ($")*CTr$”] = 0. 


(2.29) 


Summing (2.29) for j = 0,1,..., M — 1 and noticing the summation by parts formula, we have 

M — l 


/ ' \ h M ~ 1 h M_1 

E Re * (*r x - ^")Vi4($” +1 + *?)) = V E (‘T'r^r 1 - 4 E 

1=0 J E j =0 e j- o 


and 


., M — l M—l , M—l 

- y E (*" +1 )*<M**? +1 + 7 E (*")W**7 + y E W +1 )V 3 ®r 1 - 

i=o i=o j=o 

M — l M—l 

+E ^i*r l i 2 - i $ ”i 2 ) - * e ^i.i ((*r i )Vi$r i - ($”)V!$7)=o, (2.30) 

1=0 1=0 


which immediately implies (2.26). 


□ 


2-4- Error estimates 

Let 0 < T < T* with T* being the maximal existence time of the solution, and denote fix = [0,T] x 
fl. Motivated by the nonrelativistic limit of the Dirac equation [T8] and the dispersion relation (1.26), 
we assume that the exact solution of (2.1) satisfies $ £ C 3 ([0,T]; (L°°(f2)) 2 ) fl C ,2 ([0,T]; (Wd’°°(fl!)) 2 ) fl 


C 1 ([0, T\; (W 2 ’°°(n)) 2 ) n C([0,T]; (W 3 ’*^)) 2 ) and 


(^) 


Qr+s 


dt r d x s 


$ 


< -E 0<r<3, 0<r + s<3, 0 < e < 1, 


(2.31) 


L“([0,T];(L~(f2))2) 


where W’™ ,00 (fl) = {u \ u £ d l x u(a) = d x u{b), 1 = 0 ,..., m— 1} for m> 1 and here the boundary 

values are understood in the trace sense. In the subsequent discussion, we will omit fl when referring to the 
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space norm taken on SI. In addition, we assume the electromagnetic potentials V £ C(flr) and A\ £ C(SIt) 
and denote 


(B) 


Vtnax := max \V(t,x)\, A lt , 
(t,x)G^T 


:= max \A\(t, x)\. 
(t,x)EQ,T 


(2.32) 


(2.33) 


Define the grid error function e" = (eft, eft,..., eft f ) T £ Xm as: 

e" = $(t„, Xj ) - j = 0,1,. .., M, n > 0, 

with <f>” being the approximations obtained from the FDTD methods. 

For the CNFD (2.9), we can establish the following error bound (see its proof in Appendix A). 

Theorem 2.1. Under the assumptions (A) and (B), there exist constants ho > 0 and tq > 0 sufficiently 
small and independent of e, such that for any 0 < e < 1, 0 < h < ho and 0 < r < To, we have the following 
error estimate for the CNFD (2.9) with (2.10) 


^ h 2 t 2 

< -—r. 


T 

0 < n < 

T 


For the LFFD (2.6), we assume the stability condition 

e 2 h 

0 < t < 


h > 0, 0 < e < 1, 


(2.34) 


(2.35) 


£ 2 h \max + \//r 2 + e 2 (l + £AAi imax ) 2 

and establish the following error estimate (see its proof in Appendix B). 

Theorem 2.2. Under the assumptions (A) and (B), there exist constants h o > 0 and r 0 > 0 sufficiently 
small and independent of e, such that for any 0 < e < 1, when 0 < h < ho and 0 < r < To and under the 
stability condition (2.35), we have the following error estimate for the LFFD (2.6) with ( 2.1(f\) and (2.11) 

l|e n || P <- + i 0 < n < —. (2.36) 


Similarly to the proofs of the LFFD and CNFD methods, error estimates for SIFD1 (2.7) and SIFD2 


(2.8) can be derived and the details are omitted here for brevity. For the SIFD2 (2.8), we assume the 

1 


stability condition 


0 < r < 


V m 


A i, 


h > 0, 0 < e < 1, 


(2.37) 


and establish the following error estimates. 


Theorem 2.3. Under the assumptions (A) and (B), there exist constants h 0 > 0 and Tq > 0 sufficiently 
small and independent of e, such that for any 0 < e < 1, when 0 < h < ho and 0 < r < tq and under the 
stability condition (2.15), we have the following error estimate for the SIFD1 (2.1) with (2.10) and (2.11) 

■2 


U 2 




T 

0 < n < 

T 


(2.38) 


Theorem 2.4. Under the assumptions (A) and (B), there exist constants ho > 0 and r 0 > 0 sufficiently 
small and independent of s, such that for any 0 < £ < 1, when 0 < h < ho and 0 < r < To and under the 
stability condition (2.31), we have the following error estimate for the SIFD2 (2.8) with (2.10) and (2.11) 


h 2 t 2 

I P <- + -«, 

£ £ b 


T 

0 < n < —. 

T 


(2.39) 


Based on Theorems 2.1|2.4 the four FDTD methods studied here share the same temporal/spatial 
resolution capacity in the nonrelativistic limit regime. In fact, given an accuracy bound <5 > 0, the e- 
scalability of the four FDTD methods is: 


o(e 3 v^) =0(e 3 ), h = O (Vfie) = O (y/e) , 0 < £ < 1. 


(2.40) 
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3. An EWI-FP method and its analysis 


In this section, we propose an exponential wave integrator Fourier pseudospectral (EWI-FP) method to 
solve the Dirac equation (1.9) (or (1.21|) and establish its stability and convergence in the nonrelativistic 
limit regime. Again, for simplicity of notations, we shall only present the numerical method and its analysis 
for (2.1) in ID. Generalization to (1.9) and/or higher dimensions is straightforward and the results remain 
valid without modifications (see generalizations in Appendix D). 

3.1. The EWI-FP method 
Denote 

Y M = Z M x Z M , with = span|0;(a;) = e w( * _a) , / = - y,- y + 1,..., y - l| . 

Let [C p (D)] 2 be the function space consisting of all periodic vector function U(x) : O = [a, b) —> C 2 . For 
any U(x) G [C'p(D)] 2 and U € Xm, define Pm ■ [L 2 (D)] 2 —► Ym as the standard projection operator p5] , 
Im ■ [Cp(D)] 2 —> Ym and Im ■ Xm —> Ym as the standard interpolation operator f69], he. 


M/2-1 M/2-1 

(P M U)(x) = J2 ^e <w(B - o) , (I M U)(x)= J2 ^ e * w(x “ a) > a < x < b, 

l=-M/2 l—— M/2 


with 

U, = 


M-l 




0 —2ijln/M 


M M 

l = -,-h 1, 

2 2 


M 


(3.1) 


-1, (3.2) 


where Uj = U(xj) when U is a function. 

The Fourier spectral discretization for the Dirac equation (2.1) is as follows: 
Find $M(t,x) G Y M , he. 


M/2-1 

l=-M/2 


a < x < b, t > 0, 


such that for a < x < b and t > 0, 


id t $ M (t,x) = 


1 


-~<Jld x + cr 3 
£ e A 


®M{t,x) + P M {V$ M )(t,x) - <JiP M {Ai<b M )(t,x). 


(3.3) 


(3.4) 


Substituting (3.3) into (3.4), noticing the orthogonality of <pi(x), we get for l = — ..., -g— 1, 

W 




1 

(Ti + — a 3 


($ M ),(t) + (V$ M ) l (t)-<T 1 (A 1 $ M ) l (t)=0, t> 0. (3.5) 


For each l (Z = —4/,—+ 1,..., 4y — 1), when t is near t = t n {n > 0), we rewrite the above ODEs as 

d 
ds 


i -j-^(&M)i(tn + s) — (^M)i(tn + s) + Fp(s), S £ R, 


where 

Fp(s) = (G<S> M )i(tn + s), G{t,x) = V(t,x)I 2 - oiAi(t,a;) 

and F/ = meai + <73 = Qi A (Qi)* with 

l+Si _ ew 


s, t e R, 


(3.6) 

(3.7) 


- f 1 W £ \ r I V 2S ‘(i +<50 

1 - ( -1J ’ = I _-ail_ _li 


i(i+<Si) 


1+ii 


, A = 


, \/2Si(l+Si) ^25i(l+6i) 
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<5/ 0 

0 -<5; 


= a/ 1 + £2 tf- (3-8) 


























Solving the above ODE (3.6) via the integrating factor method, we obtain 

+ = e~ lsTl,e 2 (§ M )i{tn) -i [ e l{w ~ s)Tl/e 2 F{ l tw)dw 1 

Jo 

Taking s = r in (3.9) we have 

($ M )i(t n + 1 ) = e~ lTVlle2 ($ M )i(tn) - i f 


s G 


(3.9) 


1 FJ l (w)dw. (3.10) 

To obtain an explicit numerical method with second order accuracy in time, we approximate the integrals in 


/ o 
ira 

(|3.10|) via the Gautschi-type rules [40l 0HJ 09], which have been widely used for integrating highly oscillatory 

~-- r ' 3°(0), 


odes m mm saiira m, as 


i(w-T) 


r i Z?0 


Fi(w) dw 


2p— 1 


/0 


/0 


dwF?{ 0) = -i£ 2 T] 


h-e 


(3.11) 


and 


?l Ff{w)dw- 


^Fr(0) + w5^Fn0))dw 


= —*e 2 rr 1 


h~e 




m o) 


—*e 2 rr ; _1 + £ 4 r" 2 


+ £ 4 r- 2 (/ 2 -e-ft r! )]^-^"(0), n> 1 


(3.12) 


where we have approximated the time derivative dtF™(s) at s = 0 for n > 1 by the finite difference as 


WO) « «"(0) = 


F-(0) - E ; n_1 (0) 


(3.13) 


Now, we are ready to describe our scheme. Let $m( x ) be the approximation of x) (n > 0). Choosing 

®m( x ) = (Pm &o)(x), an exponential wave integrator Fourier spectral (EWI-FS) discretization for the Dirac 
equation (2.1) is to update the numerical approximation d>^ 1 (cc) G Ym (n = 0,1,...) as 


M/2-1 _ 

E Si?V w(B ~ 0) > 

l=-M/ 2 


a < x < b, n > 0, 


(3.14) 


where for Z = -f-,f - 1, 


(*2T), = 


e -irr,/ e 2 (^ ) fear- 1 


h-e 


(G(to)$ 0 M ) v 


n = 0, 


e" iT r,/ £2 (^) ; - iQ\ 1] M (G(t n )n i ) l ~ *QPV)<V {G{t n )^ M ) v n> 1, 
with G(t) denoting G(t,x) and the matrices Q^\t) and Q\ 2 \t) given as 


qPE) = ~ i£ ' 2T ~i 


2-p— 1 


In - e~^ Tl 


, q¥\t) = -ie 2 rT^ + e 4 rr 2 (/ 2 - e"S r ‘) . 


(3.15) 


(3.16) 


The above procedure is not suitable in practice due to the difficulty in computing the Fourier coefficients 
through integrals in (3.2). Here we present an efficient implementation by choosing ^^(a:) as the interpolant 
of 4>o(:r) on the grids {xj,j = 0,1,..., M } and approximate the integrals in (3.2) by a quadrature rule. 

Let be the numerical approximation of <I>(f n , x 3 ) for j = 0,1,2,...,M and n > 0, and denote 
4> n G Xm as the vector with components Choosing = $ 0 (xj) (j = 0,1,... , M), an EWI Fourier 
pseudospectral (EWI-FP) method for computing <f> n+1 for n > 0 reads 


M/ 2-1 


$ 


n+1 


= E ( $n+1 )d 


,2 ijln/M 


j = 0,1 


(3.17) 


l = — M/2 


14 
























where 


(® n+1 )i = 


e -irr,/ e 2 ^ o ^ _ i£ .2 T -l 
p -irTi/e 2 f^\. _ iC )[^ 


h-e~7* 


-5r, 


(G(io)$ 0 )p 

( 2 )/ 


n = 0, 


(* n )i - iQj i; (r)(G(i n )$n ); _ iQY>{T)8^{G{t n )^) v n > 1. 


(3.18) 


The EWI-FP (3.17)-(3.18) is explicit, and can be computed efficiently by the fast Fourier transform (FFT). 
The memory cost is O(M) and the computational cost per time step is 0[M log M). 


3.2. Linear stability analysis 

To consider the linear stability, we assume that in the Dirac equation (2.1), the external potential fields 
are constants, i.e. Ai(t,x) = A ° and V(t,x) = V° with A ° and V° being two real constants. In this case, 
we adopt the Von Neumann stability requirement that the errors grow exponentially at most. Then we have 


Lemma 3. The EWI-FP method (3.17)-(3.18) and EWI-FS method (3.1^)-(3.15) are stable under the 
stability condition 

0 < r < 1, 0 < e < 1. (3.19) 


Proof: We shall only prove the EWI-FS case (3.14)-(3.15), as the EWI-FP method case (3.18) is quite 
the same. Similarly to the proof of Lemma[l] noticing (3.151, (3.7), (3.10) and (3.12), we find that 


= Zie- irr '^(i J>°), — i e i(w-T)Ti/ e 2 (y°/ 2 _ A°ai) + 1)) 


(3.20) 


Denoting C = |V°| + |A°|, taking the l 2 norms of the vectors on both sides of (3.20) and then dividing both 
sides by the l 2 norm of ($°)j, in view of the properties of e~ lsTl ^ s , we get 


G 


G 


which implies 


Thus, we obtain 


161 - 


161 < 1 +Gt+ — T 161 + yrT, 


1 + 3Gr/2\ 2 1 + 5Gr + 9G 2 r 2 /4 (1 + 5Gr/2) 




< 


M M 

l6l<l + 4Gr, l = -- 1, 


and it follows that the EWI-FS ( |3.14 )-(3.15) is stable under the stability condition (3.19). 


(3.21) 

(3.22) 

(3.23) 

□ 


3.3. Error estimates 


In order to obtain an error estimate for the EWI methods (3.14)-(3.15) and (3.17)-(3.18), motivated by 
the results in hum], we assume that there exists an integer mo > 2 such that the exact solution &(t,x) of 
the Dirac equation (2.1) satisfies 


(G) ||^||i~([0,T];(Hj* o )2) 


< 1 , 


||dt$||z,°°([0,T] ; (Z,2)2) 


<1 


||0tt$|U«([O,T];(£2)2) 


<1 


where H k (Ll) = {zi | u £ H k (Ll), d l x u(a) = d l x u(b), l = 0 ,...,k — 1}. In addition, we assume the electro¬ 
magnetic potentials satisfy 


(-D) ||46lrV 2 ,~([0,T];Z,=°) + ll^i ||rv 2 ’ oo ([0,T];L° Ci ) ^ 1- 

The following estimate can be established (see its proof in Appendix C). 
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Theorem 3.1. Let ^^(x) be the approximation obtained from the EWI-FS (3.1f)-(3.15). Under the as¬ 
sumptions (C) and (D), there exists h 0 > 0 and tq > 0 sufficiently small and independent of e such that, 


for any 0 < e < 1, when 0 < h < ho and 0 < r < tq satisfying the stability condition (2.31), we have the 
following error estimate 


\^{t n ,x) - ^h{x )\\ L 2 < — + h m °, 0 < n < 


T 


(3.24) 


Remark 3.1. The same error estimate in Theorem 3.1 holds for the EWI-FP (3.11)-(3.18) and the proof 
is quite similar to that of Theorem\3.1\ 


4. A TSFP method and its analysis 

In this section, we present a time-splitting Fourier pseudospectral (TSFP) method to solve the Dirac 
equation( |1.9[ ) (or (1.21)) which has been proposed and studied for the Maxwell-Dirac equation [T?] |5T] . 
Again, for simplicity of notations, we shall only present the numerical method and its analysis for (2.1) in 


ID. Generalization to ( |1.9| ) and/or higher dimensions is straightforward and results remain valid without 
modifications (see generalizations in Appendix D). 


From time t = t n to time t = t n + 1 , the Dirac equation (2.1) is split into two steps. One solves first 


id t <$(t,x) = 


i o , 1 
-~criO x + 3 


<f>(f, 2 ;), x £ n, 


(4.1) 


£ £* 

with the periodic boundary condition ( |2.2| ) for the time step of length r, followed by solving 

id t <b(t,x) = [-Ai(f, 2 ;)cri + V(t,x)I 2 \ ®(t,x), x £ D, (4.2) 

for the same time step. Eq. (14.1|) will be first discretized in space by the Fourier spectral method and then 


integrated (in phase or Fourier space) in time exactly [T7J. For the ODEs (4.2), we can integrate analytically 
in time as 

4>(t, 2 ;) = e ~ l ft n l v ( s ’ x '> l2 ~ Al ( s ’ x l (Tl l ds <fr(t n ,x), a<x<b , t n <t<t n+ 1 . (4.3) 

In practical computation, from time t = t n to t = t n+ 1 , one often combines the splitting steps via the 
standard Strang splitting m - which results in a second order time-splitting Fourier pseudospectral (TSFP) 
method - as 


$ 


(i) 


M/ 2-1 


M/ 2-1 


= E 

1—-M/2 


0 -jiT ! /2e 2 


e ($") I e iw ^-“) = ® ie 


—irDi /2e 2 


(QiY (*" 


2ijl 7T 

, e M 


$( 2 ) = e -< n+l dt $( i ) = p e ~ iA 

3 2 2 


■i P* $ 

3 3 


J——M/2 
( 1 ) 


M/2-1 
1 ——M /2 


e -irV l /2e 2 ($(2))^ 




M/2-1 


j = 0, 

0 —irDi/2e 2 


n > 0, 




(4.4) 


;=—m/2 

K 1 ) _ r 4 "+1 


where f *™ +1 G(t , Xj)dt = V'j 1 ' 1 / 2 - A^j a\ = Pj A j P* with V) (1) = f *" +1 V (t, Xj)dt, A^j = A 1 (t, Xj)dt, 

A j = diag(Ay_, A J)+ ) with A y± = vj 1 ^ ± A^ 1 ], and Pj = J 2 if A^j = 0 and otherwise 


P, = P(°) : = 



(4.5) 


Remark 4.1. Again, if the definite integrals in f* n+1 A(t,Xj) dt cannot be evaluated analytically, we can 
evaluate them numerically via the Simpson’s quadrature rule as 

J A\ (t, Xj ) dt ~ ~ A-, ( t n , Xj ) H - 4Ai ^ t n -j- —, Xj ^ T A^ (fn+i 5 ) 

V(t,Xj)dt^^ V(t n ,Xj) +4V (t n -(- +V(t n+ i,Xj) 
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Lemma 4. The TSFP (4-4) conserves the mass in the discretized level, i.e. 


M -1 M — l M—l 

Il*l2> ■■= h l^l 2 = h E l^l 2 = Mp = h E l‘M^)| 2 , n > 0. (4.6) 

j=o j'=o i=o 

Proof: The proof is quite standard and similar to that of Lemma [2] We omit it here. □ 

From Lemma |4j we conclude that the TSFP (4.4) is unconditionally stable. In addition, under proper 
assumptions of the exact solution 4>(t, x) and electromagnetic potentials, it is easy to show the following 
error estimate via the formal Lie calculus introduced in l56l, 




T 

0 < n < —, 

T 


(4.7) 


where mo depends on the regularity of <f>(f, x). We omit the details here for brevity. 


5. Numerical comparison and applications 

In this section, we compare the accuracy of different numerical methods including the FDTD, EWI-FP 


and TSFP methods for the Dirac equation (1.21) in ID in terms of the mesh size h, time step r and the 


parameter 0 < e < 1. We will pay particular attention to the £-scalability of different methods in the 


nonrelativistic limit regime, i.e. 0<e<l. Then we simulate the dynamics of the Dirac equation (1.21) in 
2D with a honeycomb lattice potential by the TSFP method. 

5.1. Comparison of spatial/temporal resolution 


as 


To test the accuracy, we choose the electromagnetic potentials in the Dirac equation (1.21) with d = 1 

V(t,x) = --E ieK, t> 0, (5.8) 

1 + x z 


Ai(t,x) = 


(.x + l ) 2 


1 + x 2 


and the initial data as 


(pi{0,x) = e 


- * 2 /2 


<£ 2 ( 0 , x) = e 


- r/2 


X G 


(5.9) 


The problem is solved numerically on an interval Q = (—16,16) with periodic boundary conditions on dfl. 
The ‘reference’ solution 4>(f, x) = x), 02(L x)) T is obtained numerically by using the TSFP method 

with a small time step and a very fine mesh size, e.g. r e = 10 ' and h e = 1/16 or h e = 1/4096 for 
the comparison of the EWI-FP/TSFP methods or the FDTD methods, respectively. Denote <!>/ r as the 
numerical solution obtained by a numerical method with mesh size h and time step r. In order to quantify 
the convergence, we introduce 


eh, T {tn) = ||$" - $(W 


P = 


\ 


M—l 




3=0 


Table 5.1 lists spatial errors eh, Te (t = 2) with different h (upper part) and temporal errors e/, e)T (t = 


2) with different r (lower part) for the LFFD method (2.6). Tables [5~~2| 5.6 show similar results for the 


SIFD1 method (2.7), SIFD2 method (2.8), CNFD method (2.9), EWI-FP method (3.17)-(3.18) and TSFP 


method (4.4), respectively. For the LFFD and SIFD1 methods, due to the stability condition and accuracy 

2 eo/2 J < e < 1, 


requirement, we take 


8j{e) = 


el/V 0 <£< £ o / 2 j ', 


j = 0,1,... 


in Tables 5.1 and 5.2 For comparison, Table 5.7 depicts temporal errors of different numerical methods 
when £ = 1 for different r, Table [5~8] depicts temporal errors of different numerical methods under different 
£-scalability. 
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Table 5.1: Spatial and temporal error analysis of the LFFD method for the Dirac equation (1.21) in ID. 


Spatial Errors 

00 

i—1 

II 

O 

ho/2 

h 0 /2 2 

h 0 /2 3 

ho/ 2 4 

£o = 1 

1.06E-1 

2.65E-2 

6.58E-3 

1.64E-3 

4.10E-4 

order 

- 

2.00 

2.01 

2.00 

2.00 

£o/2 

9.06E-2 

2.26E-2 

5.64E-3 

1.41E-3 

3.51E-4 

order 

- 

2.00 

2.00 

2.00 

2.00 

£o/2 2 

8.03E-2 

2.02E-2 

5.04E-3 

1.25E-3 

3.05E-4 

order 

- 

1.99 

2.00 

2.01 

2.02 

W2 3 

9.89E-2 

2.47E-2 

6.17E-3 

1.54E-3 

3.85E-4 

order 

- 

2.00 

2.00 

2.00 

2.00 

eo/2 3 

9.87E-2 

2.48E-2 

6.18E-3 

1.54E-3 

3.83E-4 

order 

- 

2.00 

2.00 

2.00 

2.01 

Temporal Errors 

TO = 0.1 
ho = 1/8 

to/8 

ho/8Si (e) 

to/8 2 

h 0 /8 2 5 2 (e) 

t 0 /8 3 

ho/8 3 6 3 (e) 

to/8 4 

h 0 /8 4 5 i {e) 

£o = 1 

1.38E-1 

1.99E-3 

3.11E-5 

4.86E-7 

7.59E-9 

order 

- 

2.04 

2.00 

2.00 

2.00 

£ o/2 

unstable 

1.14E-2 

1.77E-4 

2.77E-6 

4.32E-8 

order 

- 

- 

2.00 

2.00 

2.00 

eo/2 2 

unstable 

4.59E-1 

7.01E-3 

1.05E-4 

1.64E-6 

order 

- 

- 

2.01 

2.02 

2.00 

£o/2 3 

unstable 

unstable 

4.14E-1 

6.42E-3 

1.00E-4 

order 

- 

- 

- 

2.00 

2.00 

£o/2 4 

unstable 

unstable 

unstable 

4.04E-1 

6.00E-3 

order 

- 

- 

- 

- 

2.02 


From Tables |5.1||5.8[ and additional numerical results not shown here for brevity, we can draw the 
following conclusions for the Dirac equation by using different numerical methods: 


(i) . For the discretization error in space, for any fixed e = £o > 0, the FDTD methods are second-order 
accurate, and resp., the EWI-FP and TSFP methods are spectrally accurate (cf. each row in the upper 
parts of Tables |5.l| 5 lljand Table |5?7| ) . For 0 < e < 1, the errors are independent of e for the EWI-FP and 
TSFP methods (cf. each column in the upper parts of Tables 5.5|5.6 ), and resp., are almost independent of 
e for the FDTD methods (cf. each column in the upper parts of Tables 5.1|5.4 ). In general, for any fixed 
0 < e < 1 and h > 0, the EWI-FP and TSFP methods perform much better than the FDTD methods in 
spatial discretization. 

(ii) . For the discretization error in time, in the 0(1) speed-of-light regime, i.e. £ = 0(1), all the 
numerical methods including FDTD, EWI-FP and TSFP are second-order accurate (cf. the first row in the 
lower parts of Tables 5.1||5.6 ). In general, the EWI-FP and TSFP methods perform much better than the 
FDTD methods in temporal discretizations for a fixed time step. In the non-relativistic limit regime, i.e. 
0 < £ <C 1, for the FDTD methods, the ‘correct’ £-scalability is r = 0 (e 3 ) which verifies our theoretical 
results; for the EWI-FP and TSFP methods, the ‘correct’ £-scalability is r = 0 (e 2 ) which again confirms 
our theoretical results. In fact, for 0 < £ < 1, one can observe clearly second-order convergence in time for 
the FDTD methods only when r < £ 3 (cf. upper triangles in the lower parts of Tables 5.l|5.4 ), a nd resp. , for 
the EWI-FP and TSFP methods when r < £ 2 (cf. upper triangles in the lower parts of Tables 5.5||5.6 ). In 
general, for any fixed 0 < £ < 1 and r > 0, the TSFP method performs the best, and the EWI-FP method 


performs much better than the FDTD methods in temporal discretization (cf. Table 5.8 1 . 

(iii). From Table 5.6 our numerical results suggest the following error bound for the TSFP method when 
r < £ 2 , 

||S(im •) - /m($") iu= < h m ° + 0 < n < -, 

r 
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(5.10) 


































Table 5.2: Spatial and temporal error analysis of the SIFD1 method for the Dirac equation (1.21) in ID. 


Spatial Errors 

00 

i—1 

II 

O 

ho/2 

ho/2 2 

ho/2 3 

ho/ 2 4 

eo = 1 

1.06E-1 

2.65E-2 

6.58E-3 

1.64E-3 

4.10E-4 

order 

- 

2.00 

2.01 

2.00 

2.00 

£o/2 

9.06E-2 

2.26E-2 

5.64E-3 

1.41E-3 

3.51E-4 

order 

- 

2.00 

2.00 

2.00 

2.00 

£o/2 2 

8.03E-2 

2.02E-2 

5.04E-3 

1.25E-3 

3.05E-4 

order 

- 

1.99 

2.00 

2.01 

2.02 

W2 3 

9.89E-2 

2.47E-2 

6.17E-3 

1.54E-3 

3.85E-4 

order 

- 

2.00 

2.00 

2.00 

2.00 

eo/2 3 

9.87E-2 

2.48E-2 

6.18E-3 

1.54E-3 

3.83E-4 

order 

- 

1.99 

2.00 

2.00 

2.01 

Temporal Errors 

TO = 0.1 

h 0 = 1/8 

to/8 

ho/ 8$i(e) 

to/8 2 

h 0 /8 2 S 2 {£) 

t 0 /8 3 

ho/8 3 S 3 (£) 

to/8 4 

/i 0 /8 4 <5 4 (e) 

£o = 1 

1.44E-1 

2.09E-3 

3.27E-5 

5.11E-7 

7.98E-9 

order 

- 

2.03 

2.00 

2.00 

2.00 

£ o/2 

unstable 

2.99E-2 

4.67E-4 

7.30E-6 

1.14E-7 

order 

- 

- 

2.00 

2.00 

2.00 

e 0 /2 2 

unstable 

8.18E-1 

1.54E-2 

2.41E-4 

3.77E-6 

order 

- 

- 

1.91 

2.00 

2.00 

£o/2 3 

unstable 

unstable 

7.99E-1 

1.31E-2 

2.05E-4 

order 

- 

- 

- 

1.98 

2.00 

£o/2 4 

unstable 

unstable 

4.19E-1 

7.97E-1 

1.26E-2 

order 

- 

- 

- 

-0.31 

1.99 


which is much better than (4.7) for the TSFP method in the nonrelativistic limit regime. Rigorous mathe¬ 
matical justification for (5.10) is on-going. 

From Tables |5.1||5.4| in the numerical example, we could not observe numerically the e-dependence in 
the spatial discretization error for the FDTD methods, i.e. 4 in front of h 2 , which was proven in Theorems 
2.1||2.4 In order to investigate the spatial e-resolution of the FDTD methods, we consider the Dirac equation 


(2.1) onH= (—1,1) with no electromagnetic potential - the free Dirac equation, i.e. 

Ai(f,:r) = 0, V(t, x) = 0, x £ (— 1,1), t > 0. 


The initial data in (2.2) is taken as 


(0, x) = e 9 ™ ( * +1) , 0 2 (O, x) = e 9 ™ ( * +1) , 


-1 < x < 1. 


(5.11) 


(5.12) 


Table 5.9 shows the spatial errors e/i iTe (t = 2) of the CNFD method with different h. The results for the 
LFFD, SIFD1 and SIFD2 methods are similar and they are omitted here for brevity. From Table [/TO] we 
can conclude that the error bounds in the Theorems |2.l|[2~f| are sharp. 


Based on the above comparison, in view of both temporal and spatial accuracies and resolution capacity, 
we conclude that the EWI-FP and TSFP methods perform much better than the FDTD methods for 
the discretization of the Dirac equation, especially in the nonrelativistic limit regime. For the reader’s 
convenience, we summarize the properties of different numerical methods in Table |5.10| 

As observed in [QB 01) , the time-splitting spectral (TSSP) method for the Schrodinger equation performs 
much better for the physical observable, e.g. density and current, than for the wave function, in the 
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Table 5.3: Spatial and temporal error analysis of the SIFD2 method for the Dirac equation (1.21) in ID. 


Spatial Errors 

00 

1 —1 

II 

O 

ho/2 

h 0 /2 2 

h 0 /2 3 

ho/2 4 

£o = 1 

1.06E-1 

2.65E-2 

6.58E-3 

1.64E-3 

4.10E-4 

order 

- 

2.00 

2.01 

2.00 

2.00 

£o/2 

9.06E-2 

2.26E-2 

5.64E-3 

1.41E-3 

3.51E-4 

order 

- 

2.00 

2.00 

2.00 

2.00 

£o/2 2 

8.03E-2 

2.02E-2 

5.04E-3 

1.25E-3 

3.05E-4 

order 

- 

1.99 

2.00 

2.01 

2.02 

W2 3 

9.89E-2 

2.47E-2 

6.17E-3 

1.54E-3 

3.85E-4 

order 

- 

2.00 

2.00 

2.00 

2.00 

eoT? 

9.87E-2 

2.48E-2 

6.18E-3 

1.54E-3 

3.83E-4 

order 

- 

1.99 

2.00 

2.00 

2.01 

Temporal Errors 

To=0.1 

to/8 

(N 

OO 

p 

ro/8 3 

OO 

1? 

£o = l 

1.72E-1 

2.59E-3 

4.05E-5 

6.33E-7 

9.89E-9 

order 

- 

2.01 

2.00 

2.00 

2.00 

£o/2 

1.69 

3.57E-2 

5.58E-4 

8.72E-6 

1.36E-7 

order 

- 

1.86 

2.00 

2.00 

2.00 

W2 2 

2.59 

8.66E-1 

1.63E-2 

2.55E-4 

3.98E-6 

order 

- 

0.52 

1.91 

2.00 

2.00 

£o/2 3 

2.67 

2.89 

8.43E-1 

1.37E-2 

2.14E-4 

order 

- 

-0.04 

0.59 

1.98 

2.00 

e o/2 4 

3.07 

3.56 

5.19E-1 

8.37E-1 

1.28E-2 

order 

- 

-0.07 

0.93 

-0.23 

2.01 


semiclassical limit regime with respect to the scaled Planck constnat 0<e<l. In order to see whether this 
is still valid for the TSFP method for the Dirac equation in the nonrelativistic limit regime, let p n = r | 2 , 

J ra = 1(<I>™ T with T the numerical solution obtained by the TSFP method with mesh size h and 

time step t, and define the errors 


e h p ’ T (tn) := II p n -p(t ri 


N— 1 

I * 1 = h Yl I Pj-Ptfnt- 

J—0 


e h j' T (t n ) := \\3 n -J(t n ,-)\\p=hy]\3]-J(t n ,x j )\. 


N-l 


1=0 


Table 

From 


5.11 lists temporal errors e 
this Table 


r (t = 2) and ej ,T (t = 2) with different r for the TSFP method (4.4). 


we can see that the approximations of the density and current are at the same order as 
for the wave function by us ing th e T SFP method. The reason that we can speculate is that p = 0(1) and 
J = 0(e _1 ) (see details in ( |l.ll| ) or ( |l.23[ )) in the Dirac equation, where in the Schrodinger equation both 
density and current are all at 0(1), when 0 < £ <1. 


5.2. Dynamics of the Dirac equation in 2D 

Here we study numerically the dynamics of the Dirac equation (1.21) in 2D with a honeycomb lattice 
potential, i.e. we take d = 2 and 


x) = A 2 (t,x) = 0, H(f,x) = cos ( ^= e i ' + cos ^^= e 2 • x 


47T 

■ cos | —^=e 3 • x 


with 


e 1 = (-l,0) T , e 2 = (1/2, Vz/2) 1 


e 3 


= (1/2,-V3/2) t . 


(5.13) 

(5.14) 


The initial data in (1.22) is taken as 


2 +V± , , n s _ (^-lp + y 2 

2 i <^ 2 ( 0 , x) = e 2 


<M0,x) = e 
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x = (x, y) T G K 2 . 


(5.15) 






























Table 5.4: Spatial and temporal error analysis of the CNFD method for the Dirac equation (1.21) in ID. 


Spatial Errors 

00 

T—1 

II 

O 

*CS 

ho/2 

ho/2 2 

ho/2 3 

ho/2 3 

£o = 1 

1.06E-1 

2.65E-2 

6.58E-3 

1.64E-3 

4.10E-4 

order 

- 

2.00 

2.01 

2.00 

2.00 

£o/2 

9.06E-2 

2.26E-2 

5.64E-3 

1.41E-3 

3.51E-4 

order 

- 

2.00 

2.00 

2.00 

2.00 

£o/2 2 

8.03E-2 

2.02E-2 

5.04E-3 

1.25E-3 

3.05E-4 

order 

- 

1.99 

2.00 

2.01 

2.02 

W2 3 

9.89E-2 

2.47E-2 

6.17E-3 

1.54E-3 

3.85E-4 

order 

- 

2.00 

2.00 

2.00 

2.00 

eoT? 

9.87E-2 

2.48E-2 

6.18E-3 

1.54E-3 

3.83E-4 

order 

- 

1.99 

2.00 

2.00 

2.01 

Temporal Errors 

76=0.1 

to/8 

to/8 2 

t 0 /8 3 

OO 

1? 

£o = 1 

5.48E-2 

8.56E-4 

1.34E-5 

2.09E-7 

3.27E-9 

order 

- 

2.00 

2.00 

2.00 

2.00 

£o/2 

3.90E-1 

6.63E-3 

1.77E-4 

2.77E-6 

4.32E-8 

order 

- 

1.96 

1.74 

2.00 

2.00 

W2 2 

1.79 

2.27E-1 

3.55E-3 

1.56E-5 

2.44E-7 

order 

- 

0.99 

2.00 

2.61 

2.00 

£o/2 3 

3.10 

4.69E-1 

2.06E-1 

3.22E-3 

5.03E-5 

order 

- 

0.91 

0.40 

2.00 

2.00 

e o/2 4 

2.34 

1.83 

8.05E-1 

2.04E-1 

3.19E-3 

order 

- 

0.12 

0.39 

0.66 

2.00 


The problem is solved nu meri call y on f1 = [—10,10] 2 by the TSFP method with mesh size h = 1/16 and 
time step r = 0.01. Figures 5.1 and 5.2 depict the densities pj(t,x) = \(f>j(t,x)\ 2 (j = 1,2) for e = 1 and 
e = 0.2, respectively. 

From Figures |5.1||5.2[ we find that the dynamics of the Dirac equation depends significantly on e. In 
addition, the TSFP method can capture the dynamics very accurately and efficiently. 


6. Conclusion 

Three types of numerical methods based on different time integrations were analyzed rigorously and 
compared numerically for simulating the Dirac equation in the nonrelativistic limit regime, i.e. 0 < e C 1 
or the speed of light goes to infinity. The first class consists of the second order standard FDTD methods, 
including energy conservative/ nonconservative and implicit/semi-implicit/explicit ones. In the nonrelativis¬ 
tic limit regime, the error estimates of the FDTD methods were rigorously analyzed, which suggest that 
the e-scalability of the FDTD methods is r = 0(e 3 ) and h = 0(y/e). The second class applies the Fourier 
spectral discretization in space and Gautschi-type integration in time, resulting in an EWI-FP method. Rig¬ 
orous error bounds for the EWI-FP method were derived, which show that the e-scalability of the EWI-FP 
method is r = 0(e 2 ) and h = 0(1). The last class combines the Fourier spectral discretization in space 
and splitting technique in time, which leads to a TSFP method. Based on the rigorous error analysis, the 
e-scalability of the TSFP method is r = 0(e 2 ) and h = 0(1), which is similar to the EWI-FP method. From 
the error analysis and numerical results, the EWI-FP and TSFP methods perform much better than the 
FDTD methods, especially in the nonrelativistic limit regime. Extensive numerical results indicate that the 
TSFP method is superior than the EWI-FP in terms of accuracy and efficiency, and thus the TSFP method 
is favorable for solving the Dirac equation directly, especially in the nonrelativistic limit regime. Finally, we 
studied the dynamics of the Dirac equation in 2D with a honeycomb lattice potential and observed some 
interesting dynamics for different e. 
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T=0 


T=0 



Figure 5.1: Dynamics of the densities pi(t,x) = x)| 2 (left) and p2(t,x.) = |</> 2 (i, x)| 2 (right) of the Dirac 
equation in 2D with a honeycomb lattice potential when e = 1. 


22 































5 


T=0 


T=0 



Figure 5.2: Dynamics of the densities pi(t,x) = x)| 2 (left) and x) = \(/>2(t, x)| 2 (right) of the Dirac 
equation in 2D with a honeycomb potential when e = 0.2. 
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Table 5.5: Spatial and temporal error analysis of the EWI-FP method for the Dirac equation (1.211 in ID. 


Spatial Errors 

O 

II 

to 

ho /2 

ho/2 2 

h 0 /2 3 

ho/ 2 4 

£o = l 

1.10 

2.43E-1 

2.99E-3 

2.79E-6 

1.00E-8 

order 

- 

2.13 

9.02 

32.74 

16.70 

£o/2 

1.06 

1.46E-1 

1.34E-3 

9.61E-7 

5.90E-9 

order 

- 

2.69 

10.44 

37.34 

12.76 

£o/2 2 

1.11 

1.43E-1 

9.40E-4 

5.10E-7 

7.02E-9 

order 

- 

2.79 

12.33 

42.93 

8.52 

W2 3 

1.15 

1.44E-1 

7.89E-4 

3.62E-7 

6.86E-9 

order 

- 

2.83 

13.51 

46.69 

7.26 

W2 3 

1.18 

1.45E-1 

7.63E-4 

2.91E-7 

8.46E-9 

order 

- 

2.85 

13.79 

51.21 

5.86 

Temporal Errors 

7-0=0.1 

to/4 

IO 

t-o/4 3 

t-o/4 4 

£o = l 

1.40E-1 

8.51E-3 

5.33E-4 

3.34E-5 

2.09E-6 

order 

- 

2.02 

2.00 

2.00 

2.00 

£o/2 

4.11E-1 

2.37E-2 

1.49E-3 

9.29E-5 

5.81E-6 

order 

- 

2.06 

2.00 

2.00 

2.00 

W2 2 

6.03 

1.88E-1 

1.18E-2 

7.38E-4 

4.62E-5 

order 

- 

2.50 

2.00 

2.00 

2.00 

£o/2 3 

2.21 

3.98 

1.60E-1 

1.01E-2 

6.31E-4 

order 

- 

-0.42 

2.32 

2.00 

2.00 

e o/2 4 

2.16 

2.09 

3.58 

1.53E-1 

9.69E-3 

order 

- 

0.02 

-0.39 

2.27 

1.99 


Appendix A. Proof of Theorem |2.1| for the CNFD method 
Define the local truncation error = (£q,£", • ■ • ,£m) T € of the CNFD (2.9) with (2.10) as 

t-n „ -I I * S x ^(t n + 1 , Xj) + 5 x ^?(t n , Xj) 1 ^(tn+l, Xj) + &{t n , Xj) 

S.7 -= ld t + -a 1- - - ~^a 3 -2- 

+ [^i(*n+i/ 2 ,!Ej)q-i ~ b r (t n+ i/ 2 ,a.- i )J 2 ] + $(*n,a-j) , 0 < j < M - 1, n > 0. (A.l) 


Applying the Taylor expansion in (A.l), noticing (2.1) and the assumptions (A) and (B), and using the 
triangle inequality, for 0 < e < 1, we obtain 

7-2 ^_2 7-2 / [ \ 

I €j I — + y (^72 + Knax + ^4i, max J II^IICf2 r ) 

j = 0,1,..., M — 1, n > 0, (A.2) 


< 


r 2 r 2 7,2 

-\ - j < + 

£ £° £ 4 £° £ 


which immediately implies 

irn^= 0< max_ i |e”|<y + y, \\Cy<\\C\\l~<^ + Y> 0<s<l. (Ay 

Subtracting (2.9) from (A.l), noticing ( 2.33| ), we get for n > 0 


ifit e? = --a!4e” +1/2 


i 

with e 


+ y ^ e ; +1/2 + (V / +1/2 / 2 - A"j 1/2 aij e ” +1/2 + £?, 0 < j < Af - 1, (A.4) 


+1/2 _ e " +1 + e j 


for j = 0,1,..., M. and the boundary and initial conditions are given as 


““ c Mi 


--1 — c M-1j 


n> 0, e° = 0, J = 0,1,...,M. 
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Table 5.6: Spatial and temporal error analysis of the TSFP method for the Dirac equation (1.21) in ID. 


Spatial Errors 

CM 

II 

o 


ho/2 

h 0 /2 2 

ho/2 3 


ho/2 4 

eo = 1 

1.10 


2.43E-1 

2.99E-3 

2.79E-6 


9.45E-9 

order 

- 


2.13 

9.01 

32.74 


17.18 

£o/2 

1.06 


1.46E-1 

1.34E-3 

9.61E-7 


5.57E-9 

order 

- 


2.69 

10.44 

37.34 


13.14 

£o/2 2 

1.11 


1.43E-1 

9.40E-4 

5.10E-7 


6.50E-9 

order 

- 


2.79 

12.33 

42.93 


8.86 

W2 3 

1.15 


1.44E-1 

7.89E-4 

3.62E-7 


6.84E-9 

order 

- 


2.83 

13.51 

46.69 


7.27 

eo/2 3 

1.18 


1.45E-1 

7.62E-4 

2.88E-7 


7.49E-9 

order 

- 


2.85 

13.79 

51.44 


6.20 

W2 5 

1.19 


1.46E-1 

7.53E-4 

2.59E-7 


7.96E-9 

order 

- 


2.85 

13.92 

53.92 


5.70 

W2 5 

1.20 


1.47E-1 

7.49E-4 

2.63E-7 


6.90E-9 

order 

- 


2.86 

14.01 

53.37 


6.17 

Temporal Errors 

<=? 

II 

o 

t 0 /4 

ro/4 2 

r o/4 3 

t 0 /4 4 t 0 /4 5 

to/4 6 

so = 1 

2.17E-1 

1.32E-2 

8.22E-4 

5.13E-5 

3.21E-6 2.01E-7 

1.26E-8 

order 


2.02 

2.00 

2.00 

2.00 2.00 

2.00 

£o/2 

1.32 

6.60E-2 

4.07E-3 

2.54E-4 

1.59E-5 9.92E-7 

6.20E-8 

order 

- 

2.16 

2.00 

2.00 

2.00 2.00 

2.00 

e o/2 2 

2.50 

3.33E-1 

1.68E-2 

1.04E-3 

6.49E-5 4.06E-6 

2.54E-7 

order 

- 

1.45 

2.15 

2.00 

2.00 2.00 

2.00 

eo/2 3 

1.79 

1.97 

8.15E-2 

4.15E-3 

2.57E-4 1.60E-5 

1.00E-6 

order 

- 

-0.07 

2.30 

2.14 

2.01 2.00 

2.00 

£o/2 4 

1.35 

8.27E-1 

8.85E-1 

2.01E-2 

1.03E-3 6.35E-5 

3.97E-6 

order 

- 

0.35 

-0.05 

2.73 

2.14 2.01 

2.00 

W2 5 

8.73E-1 

2.25E-1 

2.33E-1 

2.49E-1 

4.98E-3 2.55E-4 

1.58E-5 

order 

- 

0.98 

-0.03 

-0.05 

2.82 2.14 

2.01 

Table 5.7: Comparison of temporal errors of different methods for the Dirac equation ( 

1.21 

) with e = 1 . 

£ = 1 

T 0 =0.1 

to/4 

ro/4 2 

t 0 /4 3 

tq/4 4 


t 0 /4 5 

LFFD 

1.38E-1 

8.00E-3 

4.98E-4 

3.11E-5 1.94E-6 

1.21E-7 

order 

- 

2.05 

2.00 

2.00 

2.00 


2.00 

SIFD1 

1.44E-1 

8.85E-3 

5.53E-4 

3.27E-5 2.16E-6 

1.35E-7 

order 

- 

2.01 

2.00 

2.04 

1.96 


2.00 

SIFD2 

1.72E-1 

1.17E-2 

7.30E-4 

4.05E-5 2.85E-6 

1.78E-7 

order 

- 

1.94 

2.00 

2.09 

1.91 


2.00 

CNFD 

5.48E-2 

3.49E-3 

2.18E-4 

1.34E-5 8.38E-7 

5.23E-8 

order 

- 

1.99 

2.00 

2.01 

2.00 


2.00 

EWI-FP 

1.40E-1 

8.51E-3 

5.33E-4 

3.34E-5 2.09E-6 

1.30E-7 

order 

- 

2.02 

2.00 

2.00 

2.00 


2.00 

TSFP 

1.32E-2 

8.22E-4 

5.13E-5 

3.21E-6 2.01E-7 

1.26E-8 

order 

- 

2.00 

2.00 

2.00 

2.00 


2.00 


Similarly to the proof for Lemma, multiplying (A.4) from the left by hr ( e j +1//2 ) 5 taking the imaginary 
part, then summing for j = 0,1,..., M — 1, using the triangle inequality and Young’s inequality, noticing 
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Table 5.8: Comparison of temporal errors of different numerical methods for the Dirac equation (1.21) under 
proper £-scalability. 


r = 0(e 3 ) 
t = 0(h) 

£0 = 1 

£o/2 

£ o/2 2 

£ o/2 3 

£o/2 4 

ho = 1/8 

ho/2 

ho/2 2 

ho/2 3 

ho/2 4 

To = 0.1 

t 0 /8 

to/8 2 

to/8 3 

to/8 4 

LFFD 

1.38E-1 

1.14E-2 

7.01E-3 

6.42E-3 

6.00E-3 

SIFD1 

1.44E-1 

2.99E-2 

1.54E-2 

1.31E-2 

1.26E-2 

r = 0(e 3 ) 

£o = 1 

£ o/2 

£ o/2 2 

e o/2 3 

e o/2 4 

T 0 = 0.1 

to/ 8 

to/8 2 

to/8 3 

to/8 4 

SIFD2 

1.72E-1 

3.57E-2 

1.63E-2 

1.37E-2 

1.28E-2 

CNFD 

5.48E-2 

6.63E-3 

3.55E-3 

3.22E-3 

3.19E-3 

r = 0(e 2 ) 

£o = 1 

£o/2 

£ o/2 2 

£ o/2 3 

eo/2 4 

TO =0.1 

to/4 

t 0 /4 2 

t 0 /4 3 

to/4 4 

EWI-FP 

1.40E-1 

2.37E-2 

1.18E-2 

1.01E-2 

9.69E-3 

TSFP 

1.32E-2 

4.07E-3 

1.04E-3 

2.57E-4 

6.35E-5 

Table 5.9: 

Spatial error analysis of the CNFD method for the free Dirac equation with different h. 

£ 

£o — 1 

£ o/2 

£ o/2 2 

e o/2 3 

e o/2 4 

h 0 = 1/256 

1.61E-1 

3.21E-1 

6.35E-1 

1.21 

2.07 

ho/2 

4.03E-2 

8.05E-2 

1.59E-1 

3.07E-1 

5.43E-1 

h 0 /2 2 

1.01E-2 

2.01E-2 

3.99E-2 

7.69E-2 

1.36E-1 

h 0 /2 3 

2.52E-3 

5.03E-3 

9.97E-3 

1.92E-2 

3.41E-2 

h 0 / 2 4 

6.30E-4 

1.26E-3 

2.47E-3 

4.95E-3 

8.64E-3 


Table 5.10: Comparison of properties of different numerical methods for solving the Dirac equation with M 
being the number of grid points in space. 


Method 

LFFD 

SIFD1 

SIFD2 

CNFD 

EWI-FP 

TSFP 

Time symmetric 

Yes 

Yes 

Yes 

Yes 

No 

Yes 

Mass conservation 

No 

No 

No 

Yes 

No 

Yes 

Energy conservation 

No 

No 

No 

Yes 

No 

No 

Dispersion Relation 

No 

No 

No 

No 

No 

Yes 

Unconditionally stable 

No 

No 

No 

Yes 

No 

Yes 

Explicit scheme 

Yes 

No 

No 

No 

Yes 

Yes 

Temporal accuracy 

2nd 

2nd 

2nd 

2nd 

2nd 

2nd 

Spatial accuracy 

2nd 

2nd 

2nd 

2nd 

Spectral 

Spectral 

Memory cost 

O(M) 

O(M) 

0(M ) 

O(M) 

O(M) 

O(M) 

Computational cost 

O(M) 

O(M) 

0{M In M) 

» O(M) 

0{M In M) 

0(M In M) 

Resolution 

h = 0(y/e) 

h = 0 (^/e) 

h = 0(y/e) 

h = 0(y/e) 

h = 0(1) 

h = 0(1) 

when 0<£< 1 

t = 0(e 3 ) 

T = 0(£ 3 ) 

T = 0(£ 3 ) 

r = 0(£ 3 ) 

r = 0(£ 2 ) 

r = 0(£ 2 ) 


(1.3), (A.3) and (A.5), we get 


M—l 


|e" +1 ||f 2 -|!e"||f 2 < 


m (K +1 | + |e J "|)<r(mf 2 + 


3=0 

< r(||e" +1 ||f 2 


+ ll eIl ll? 2 ) + T [ — + 


* n+1 \\l 


n > 0. 


I|e n ||?0 


(A.6) 
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Table 5.11: Temporal errors for density and current of the TSFP for the Dirac equation (1.211 in ID. 


e h D ’ T (t = 2) 

II 

O 

To/ 4 

to/4 2 

t 0 / 4 3 

tq/ 4 4 

r o/4 5 

to/ 4 6 

£0 = 1 

2.50E-1 

1.54E-2 

9.61E-4 

6.01E-5 

3.75E-6 

2.34E-7 

1.43E-8 

order 


2.01 

2.00 

2.00 

2.00 

2.00 

2.01 

£o/2 

1.22 

5.27E-2 

3.21E-3 

2.01E-4 

1.25E-5 

7.84E-7 

4.92E-8 

order 

- 

2.27 

2.02 

2.00 

2.00 

2.00 

2.01 

£ 0 / 2 ^ 

1.75 

1.86E-1 

1.00E-2 

6.20E-4 

3.87E-5 

2.42E-6 

1.52E-7 

order 

- 

1.62 

2.11 

2.01 

2.00 

2.00 

2.00 

eo/2 a 

1.11 

1.39 

2.95E-2 

1.53E-3 

9.47E-5 

5.92E-6 

3.72E-7 

order 

- 

-0.16 

2.78 

2.13 

2.01 

2.00 

2.00 

eo/2 4 

1.58 

7.58E-1 

7.81E-1 

5.46E-3 

3.01E-4 

1.87E-5 

1.17E-6 

order 

- 

0.53 

-0.02 

3.58 

2.09 

2.00 

2.00 


9.59E-1 

1.96E-1 

2.29E-1 

2.33E-1 

1.20E-3 

6.76E-5 

4.21E-6 

order 

- 

1.15 

-0.11 

-0.01 

3.8 

2.07 

2.00 

e$’ T (t = 2) 

t o =0.4 

r 0 / 4 

r 0 /4 2 

t o/4 3 

to/ 4 4 

r 0 /4 5 

r 0 /4 6 

£o = l 

1.70E-1 

1.09E-2 

6.83E-4 

4.27E-5 

2.67E-6 

1.67E-7 

1.02E-8 

order 

- 

1.98 

2.00 

2.00 

2.00 

2.00 

2.01 

e o/2 

9.15E-1 

6.39E-2 

4.00E-3 

2.50E-4 

1.56E-5 

9.76E-7 

6.08E-8 

order 

- 

1.92 

2.00 

2.00 

2.00 

2.00 

2.00 

eo/2 2 

1.58 

3.45E-1 

1.69E-2 

1.04E-3 

6.50E-5 

4.06E-6 

2.54E-7 

order 

- 

1.10 

2.18 

2.01 

2.00 

2.00 

2.00 

eo/2 3 

1.06 

1.26 

5.83E-2 

2.87E-3 

1.76E-4 

1.11E-5 

6.94E-7 

order 

- 

-0.12 

2.22 

2.17 

2.01 

2.00 

2.00 

e o/2 4 

1.11 

9.78E-1 

1.05 

2.28E-2 

1.18E-3 

7.33E-5 

4.58E-6 

order 

- 

0.09 

-0.05 

2.76 

2.13 

2.00 

2.00 

eo/2 5 

4.98E-1 

1.55E-1 

2.22E-1 

2.39E-1 

4.04E-3 

2.09E-4 

1.29E-5 

order 

- 

0.84 

-0.30 

-0.05 

2.94 

2.13 

2.01 


Summing the above inequality for n = 0,1,..., m — 1, we get 


Ie m l|? 2 -||e°||? 2 <r^| 
k—0 


e k \\f 2 +Tm 


£ £ 


2 

6 / ’ 


Taking r 0 sufficiently small, when 0 < r < r 0 , we have 


T 

0 < TO < —. 

T 


(A.7) 




^ k \\l 


rm 


h“ i \ v—> 


£ £ L 




Ti V + ^ 


k —0 x ' k -0 

Using the discrete Gronwall’s inequality and noticing ||e °||;2 = 0, we obtain 


T 

0 < TO < —. 

T 


<( h - + r " 2 


£ 


£ £ u 


T 

0 < m < —, 
r 


which immediately implies the error bound (2.34). 


(A.8) 


(A.9) 

□ 


Appendix B. Proof of Theorem |2.2| for the LFFD method 
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Define the local truncation error = (£q , f", • • •, £m) T €= Im of the LFFD (2.61 with (2.101 and (2.111 
as follows, for 0 < j < M — 1, 

~ 1 1 

ar,) + a^) - ^-cj 3 $(t„, Xj) + [A^ai - V, n / 2 ] $(t n , xj), n> 1, (B.l) 


$(0, Zj) + - ^cr 3 + V^J 2 - A^oq $o(®i)- 


1 


(B.2) 


Applying the Taylor expansion in (B.l I and (B.21, noticing (2.1) and the assumptions (A) and (B), similarly 

h 2 


to the proof of Theorem 2.1 we obtain 


1 ~ £ 4 + C 


I^I^AT + -> j = 0,1,..., M — 1, n> 1, 
J £° £ 


(B.3) 


which immediately implies 


||r || ^= o <max_i ||r|| i2 <ri|^<^ + - ( n>l, 0< £ <1. (B.4) 


Subtracting (2.6) from (B.l), noticing ( |2.33| ), we get 

iS t ej = - % ~<x 1 S x ej + ^o- 3 e" + (V^h - A”jOi) e" +1", 0 < j < M - 1, n > 1, 
where the boundary and initial conditions are given as 


-0 — 

For the first step, we have 

Denote £ n+l for n = 0,1,... as 

,n+l I 


e" = e^f, e" 3 = ej f _ 1 , n > 0, e° = 0, j = 0,l,...,M. 


~ h^ 1 T‘^ 

|e 1 || i2 =r||C°|| i2 <- I + —<- + yv. 


M— 1 


i=o 


£ n+r = | K+ 1||2 + ||e ra || 2 2 + 2 Re ( rh £ (e^V^e? ) - 21m ( ^ £ ( e ; +1 )V 3 e” 


M —1 


(B.5) 

(B.6) 

(B.7) 

(B.8) 


i=o 


and under the stability condition (2.35), e.g., r < 


e 2 rih 


f < t and 75 < |, using Cauchy inequality, we can get that 


e 2 ftV m ax + •y/h 2 +e 2 (l+e/i J 4i im ax) s 


with n = t, which implies 


2 6 


n+l||2 


It follows from (B.7) that 


l + ||e"||? a ) < £ n+1 < - (||e n+1 ||fa + ||e”|| 2 2 ) , n > 0. 


f 1 < — + ^ 


(B.9) 


(B.10) 


Similarly to the proof of Theorem 2.1 multiplying (B.5) from the left by 2 hr (e" +1 + e” - 1 ) , taking the 
imaginary part, then summing for j = 0,1,..., M — 1, using Cauchy inequality, noting ( B.4| ) and (B.9), we 
get for n > 1, 


M-l 


£” +1 - <hr Y, (Wmax + ^max) |e^ | + |Q*|) (\e] +1 \ + le^ 1 )) 

3 = 0 

/ h 2 r 2 \ 2 

<r(£"+ 1 +n+T^- + ^J , «>0. 
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Summing the above inequality for n = 1, 2,..., m — 1, we get 


£ m -£ 1 < 


mr 


k =1 


h 2 


T 

1 < m < —. 

T 


(B.ll) 


Taking To sufficiently small, using the discrete Gronwall’s inequality and noticing (B.lOl, we obtain from 
the above equation that 


£ m ( I r ) , 1 < m < - 


which immediately implies the error bound (2.36) in view of (B.9). 


(B.12) 

□ 


Appendix C. Proof of Theorem |3.1| for the EWI-FP method 
Define the error function e n (x) for n = 0,1,... as 

M/2-1 


( n l 1\ ' 

e l(l ■■= P M <f>(t n ,x) - $ n M {x) = Y e?e iMx ~ a) > 

eAX) ) l=- M/ 2 

Using the triangular inequality and standard interpolation result, we get 


a < x < b. 


(C.l) 


T 


\\$(t n ,x) - <&'m{x)\\ L 2 < ||$(t„,x) - P M $(tn,x) \\ L 2 + ||e"(a :)|| L 2 < h m ° + |je ra (x )|| L 2 0 < n < —, (C.2) 


which means that we only need estimate ||e” (re) || z> 2 . _ 

Define the local truncation error C(x) = E^-m /2 £[V w(x " a) £ Y M of the EWI-FP (3.15) for n > 0 


as 


2-p-i 




h ~ e 


(G(0)$(0)) ; , 
( 2 )/ 


n = 0, 


(C.3) 


[($(f n+1 )), - e ~ iTlr '/ £ ($(i„)) z +iQi 1) (r)(G(t n )$(t n )) l + iQ^ (r)8i (G(t n )<5>(t n )) v n > 1, 

where we write <F>(£) and G{t) in short for <F>(t, x) and G(t,x ), respectively. 

Firstly, we estimate the local truncation error £"(;r). Multiplying both sides of the Dirac equation (2.1) 
by and integrating over the interval (a, b), we easily recover the equations for ($(t))j, which are 

exactly the same as (3.6) with being replaced by $(f, x). Replacing with <I>(t,a;), we use the same 
notations F™(s) as in (3.7) and the time d erivat ives of ^"(s) enjoy the same properties of time derivatives 
of <I>(f,:r). Thus, the same representation (3.10) holds for (<E>(£„))/ with n > 1. From the derivation of the 
EWI method, it is clear that the error £, n (x) comes from the approximations for the integrals in (3.11) and 
(3.12), and we have 


f? = - * / (s) - i?(0))ds = -i 


t) 


U/pO 


ili-z) 


e r, 5 Sl F ; °(si) dsids, 


(C.4) 


>o Jo 


and for n > 1 


£P = - i / 


1 /* T 


o Jo 


d S2S2 F?(s 2 ) ds 2 d Sl + s / de^F?- 1 ^) dO^Q ) ds. (C.5) 

JO J 6t 


For n = 0, the above equalities imply |£,°| < |<9 Sl F ; °(si )|dsids and by the Bessel inequality and 

assumptions (C) and (D), we find 


M/2-1 


M/2-1 


U°(x)\\ 2 L 2 ={b-a) Y l £?| 2 ( b ~ a )T 2 I I Y \d Sl F?( Sl )\ 2 d Sl ds 

;=—m/2 ;=-m/2 

<r 2 f f ||9. 1 (G(s 1 )$(s 1 ))||i,2 d si ds < jj. 

Jo Jo £ 
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Similarly, for n > 1, we obtain 

M/2-1 

U n (x)\\l* = (b~a) Y |£f 

l=-M/ 2 
't-1 

pi pa pa i ^ 

<r 3 

./n -/n /n M 

L - o 


[ [ [ Y \ d s 2 s 2 Fi n is 2 )\ 2 ds 2 ds 1 ds + T 3 [ f [ s Y 1^1 SiK l (0 1 )\ 2 d6 1 d6ds 

Jo Jo Jo i _ m J 0 J 0 JOt i _ m 


< T 6 || 9 tt (G(t)$(t))|||ooQQ jT ];( i 2 ) 2 ) < ^g, 
where we have used the assumptions (C) and (D). Hence, we derive that 

H£°(z)Ul 2 < ^, ir(x)iiL 2 <i n> i. 


(C.6) 


(C.3), we obtain the equation for the error vector function as 


Now, we look at the error equations. For each fixed l = —M/2,...,M/2 — 1, subtracting (3.15) from 

(C.7) 


e? = 0, 


e/ +1 =e-^/ £2 er + i?r+£\ 1 < n < - - 1, 


M/2-1 ___ 

where R n (x) = ^ R n e im{x-a) g y M for n > 1 is given by 

1—-M/2 

R? = -*Q« (1) (r) - iQ¥\t) (« t -(G(0$(g), - <5r(G(*3Sf),) • (C.8) 

Using the properties of the matrices Q^(t) and Q| 2 ^(t), it is easy to verify that 


ll^ (1) (r)||2<T, ||Q[ 2) (r)|| 2 <y, l = !, 


(C.9) 


where ||Q ||2 denotes the l 2 norm of the matrix Q. Combining (C.9), (C.8) and the assumption (D), we get 

M/ 2-1 M/ 2-1 . ^ 

ll-R"(*)lli=i =(6-a) Y \R?\ 2 <(b~a)r 2 Y (| (*(£))< - ® *| + |(*7Ci)), - (if 1 )* 

1——M/2 1——M/2 ' 


k—n—1 


<r 2 £ ||*(f*,z) - ^m(x )\\ 2 L 2 < T 2 h 2m ° + r 2 ||e n (a ;)|| 2 2 + r 2 ||e" _1 (a;)|| 2 2 . 

0 


Multiplying both sides of (|C.7| by (e/ +1 + e * Tr ‘/ r e 
Cauchy inequality, we obtain 


(C.10) 

from left, taking the real parts and using 


|e/ +1 | 2 - |e/| 2 < r 

|e/ +1 | 2 + |e/| 2 

|t | 1 t 1 — 



\ R n\2 1^,2 


(C.ll) 


Multiplying the above inequality by b — a and summing together for l = —M/2,..., M/2 — 1, in view of the 
Bessel inequality, we obtain 


^ +1 {x)\\ 2 l2 - ||e"0r)|| 2 2 <r(||e” +1 (x)|| 2 2 + |K(x)|| 2 2 ) + ^||i?"(x)|| 2 2 + ^\\C(x)\\h, n> 1. (C.12) 
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Summing (C.12I for n = 1,...,m — 1, using (C.lOl and (C. 6 ), we derive 


|e m (aO||iH!eV)|lL;^E 

k=1 


\ ek ( x ) I 


mr 


L 2 


■ rnrh 


2 ran 


T 

1 < m < —. 

T 


(C.13) 


Since ||e°(a ;)|| i 2 = 0 and ||e 1 (x)||i 2 < ^ < ~y, the discrete Gronwall’s inequality will imply that for 
sufficiently small r, 




1 < m < —. 

T 


Combining (C.2) and (C.14), we draw the conclusion (3.24). 


(C-14) 

□ 


Appendix D. Extensions of the EWI-FS (3.14)-(3.15) and TSFP (4.4) in 2D and 3D 


The EWI-FS (3.14)-(3.15), EWI-FP (3.17)-(3.18) and TSFP (kL4| can be easily extended to 2D and 3D with 


tensor grids by modifying the matrices Tj in (3.8) and G(t,x) in ( |4.5| ) in the TSFP case. For the reader’s 
convenience, we present the modifications of T; in (3.8) and G(t,x) i n (|4.5| ) in 2D and 3D as follows. 

For the Dirac equation (1.21) in 2D , i.e. we take d = 2 in (1.21). The problem is truncated on 
D = (ai, b\) x ( 02 , 62 ) with mesh sizes hi = (hi — a{)/Mi and h 3 = (62 — 02 )/M 2 (Mi, M 2 two even positive 
integers) in the x- and y-direction, respectively. The wave function 4> is a two-component vector, and the 
matrix T; in (3.8) will be replaced by 


1 jk — 


1 


M 1] 


■ ( 2 ) 

■ 


e 4 1} 


■ (2) N 


-1 


( 1 ) 2 jtt 

Fj = T -: 

J h - ai 


(2) 2fc7T 
^ =5^02’ 


(D.l) 


where — ^ < j < ^ — 1, — ^ < k < ^ — 1, and the Schur decomposition Tjk = QjkDjkQ% i s given as 


2 

/. 


Qjk 


l~t~<5jfc 




y/ 26 jk(l+Sjk) y/ 26 jk(l+Sjk) 

e <4 1 )+ w4 2) i+Sjk 

V \/20, k ( I -rd :! k) ^25jk(l+Sjk) 


Djk — 


Sjk 

0 


0 

—djk 


S jk = Jl+e^y + e^>Y. (D.2) 


( 2 ) \2 


The matrix J) " +1 G(t,x)dt in (4.5) becomes f t " +1 G(t,x)dt and the Schur decomposition J t n+1 G(t,x.)dt = 


F X A X P* with 


< (1) = Itn +1 V(t,x)dt, = f*" +1 Ai(t,x)dt for l = 1,2, A^ = \J\A ( ^| 2 ■ 

diag(A x A Xj+ ), A x ,± = 14 (1) 
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\A 


± A x ’, and P x = I 2 if A x : ^ = 0 and otherwise 


2,xl 


P x = 


1 

V2 

AA+<- 
V V 2 


1 


( 1 ) 
2 ,x 


A x = 


(D.3) 


V2\ 
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For the Dirac equation ]1.9\ ) in 3D, i.e. we take d = 3 in (1.9). The problem is truncated on fl = 
(ai, 6 i) x (a 2 , 6 2 ) x (03,63) with mesh sizes hi = (hi - oi)/Mi, h 2 = (62 - a 2 )/M 2 and h 3 = (63 - a 3 )/M 3 
(Mi, M 2 , M 3 three even positive integers) in x-, y- and ^-direction, respectively. The wave function T is a 
four-component vector, and the matrix T; in (3.8) will be replaced by F jki as: 
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r jkl ~ 
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ey 
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~ £ hi 
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(D.4) 


where <j<^- 1 ,-^ < k < - 1 ,-^ < l < ^ - 1 and 


d‘’ = 


2j7T 

61 - Oi ' 


( 2 ) 

/4 = 


2kir 

62 - 02 ’ 
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(3) 

W = 


2lir 

63 - 03 


(D.5) 
























































The eigenvalues of Tjki are 


SjkhSjkl, Sjki-, Sjki, with Sjki \/ 1 "b £ 

The corresponding eigenvectors are 
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Then the Schur decomposition Tju = QjkiDjkiQ*jki is given as 
Dju — diag(<5jfc/, Sjki, Sjki, Sjki), Qjki — 


( 


V (D (2) (3) (4) 

V AU1 1 V V U1 1 V oT7l v 




The matrix f t " +1 G(t,x)dt in (4.5) becomes J t " +1 G(t,x)dt and the Schur decomposition f t " +1 G(t,x)dt = 
P X A X P; with U x (1) = //; +1 V(i, x)dt, Ag = j £ +1 Ai(t, x)di for 7 = 1,2,3, A« = ^/|a £>| 2 +|Agp + |A«p, 
A x = diag(A Xi _, A Xi _, A Xi+ , A Xj+ ), A Xj ± = V^ 1 ’ ± A x \ and P x = I 4 if A x ^ = 0 and otherwise 
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For the Dirac equation (1.9) in 2D, we simply let pg = 0, A 3 (f, x) = 0 in the above 3D case; and for 
the Dirac equation (1.9) in ID, we let /xg = /if' 1 = 0, A 2 (i, x) = A 3 (f,x) = 0 in the above 3D case. Then 
the EWI-FP (3.17)-(3.18) and TSFP (4.4) can be designed accordingly for the Dirac equation (1.9) in 2D 
and ID. 
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